/*
	etopo5.


	et5.bat = cl /AL etopo5.c proj1



*/
#include <stdio.h>
#include <string.h>
#include <math.h>

#define DEG_RAD 0.01745329251				/* converts degrees to radians */
#define R      6367.65
#define E_RAD 6400
#define E_CIR 40212
#define ROW 2160
#define COL 4320
#define MAX_V 30000

double Rho=R;
struct control
{
	char proj;
	int pnum;
	double lon_0,lon_1,lon_2,lon_p;
	double lat_0,lat_1,lat_2,lat_p;
	double pixsize,scale;
	double x0,y0;
	double c,C,F,H,k,k0,k1,k2,n,P,rho_0;
	double omega; /* elevation (tilt of surface plane) */
	double gamma; /* azimuth (east of north) */
	double horiz; /* radius of horizon arc */
}Map;
char *Path[10];
int NumPath=0;
int Elev[COL],E16[COL];
unsigned char E8[COL];
char CBuff[6481];
FILE *FpTxt;

int	wf(double lat1,double lat2,double lon1,double lon2,double r,
			double zscale,int samp,double thick,char *oname);
int scan_height(double lat1,double lat2,double lon1,double lon2,
				long *max,long *min);
int scan_height_sample(double lat1,double lat2,double lon1,double lon2,
				long *max,long *min,double sample);
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 auto_set_constants(struct control *cc);
int xy2ll(double x,double y,double *lat,double *lon,struct control cc);
int ll2xy(double lat,double lon,double *x,double *y,struct control cc);
int get_constants(struct control *cc,FILE *fptxt);
int choose_projection(struct control *cc,FILE *fptxt);
int sample(double lat1,double lat2,double lon1,double lon2,
		double maxz,double minz,double sample_rate,char *name);
FILE *fopen_path(char *filename,char *type);

main(int argc,char *argv[])

{
	int i,j,k;
	int row,col,samp,len;
	char ans,string[100],cr=0xd,oname[100];
	double lat1,lat2,lon1,lon2,clat,r,km_in,dlon,dlat,sample_rate;
	double lt1,ln1,lt2,ln2,thick=0.5;
	long max,min;
	double maxdelx,maxdely,xmax,xmin,ymax,ymin,lat,lon,x,y;

	NumPath=argc-1;
	for(i=0;i<NumPath;i++)
	{
		Path[i]=argv[i+1];
		len=strlen(Path[i]);
		if(Path[i][len-1]!=':'&&Path[i][len-1]!='\\')
			strcat(Path[i],"\\");
	}
	do
	{
		printf("  Make a sub image in Wave front format:\n\n");
		printf("   1 -- round earth\n");
		printf("   2 -- flat/solid\n");
		printf("   3 -- flat/skin\n");
		printf("   4 -- image set -- flat/solid\n");
		printf("   5 -- image set -- flat/skin -- projected)\n");
		printf("   6 -- make a sub set for 'trigrid.c'\n");
		printf("\n   a -- check max min of set\n");
		printf("\n   x -- exit program\n\n");
		ans=getch();
		if(ans=='1')
		{
			printf("give minimum and maximum latitude\n");
			scanf("%lf%lf",&lat1,&lat2);
			printf("give minimum and maximum longitude\n");
			scanf("%lf%lf",&lon1,&lon2);
			printf("Give vertical exaggeration.\n");
			scanf("%lf",&r);
			printf("Give sample rate  (1 - 12)\n");
			scanf("%d",&samp);
			printf("Do you want double sided figure?  (y or n)\n");
			if(getch()=='y')
			{
				printf("Give thickness of figure in inches.\n\n");
				scanf("%lf",&thick);
			}
			else
				thick=-1.0;
			printf("Give output file name.\n");
			scanf("%s",oname);
			wf(lat1,lat2,lon1,lon2,Rho,r,samp,thick,oname);
			exit(0);
		}
		if(ans=='2')
		{
			printf("give minimum and maximum latitude\n");
			scanf("%lf%lf",&lat1,&lat2);
			printf("give minimum and maximum longitude\n");
			scanf("%lf%lf",&lon1,&lon2);
			printf("Give elevation scale exageration\n");
			scanf("%lf",&r);
			printf("give central latitide\n");
			scanf("%lf",&clat);
			if(scan_height(lat1,lat2,lon1,lon2,&max,&min)>=0)
			{
				printf("min = %ld  max = %ld\n",min,max);
				km_in=9.0/((lat2-lat1)*111.7);
				wf_flat(lat1,lat2,lon1,lon2,clat,r,min,max,"globe.wvf",km_in,
					(double)0.125);
			}
		}
		if(ans=='3')
		{
			printf("give minimum and maximum latitude\n");
			scanf("%lf%lf",&lat1,&lat2);
			printf("give minimum and maximum longitude\n");
			scanf("%lf%lf",&lon1,&lon2);
			printf("Give elevation scale exageration\n");
			scanf("%lf",&r);
			printf("give central latitide\n");
			scanf("%lf",&clat);
			printf("give thickness in inches.\n");
			scanf("%lf",&thick);
			if(scan_height(lat1,lat2,lon1,lon2,&max,&min)>=0)
			{
				printf("min = %ld  max = %ld\n",min,max);
				km_in=9.0/((lat2-lat1)*111.7);
				wf_flat_skin(lat1,lat2,lon1,lon2,clat,r,min,max,"globe.wvf",km_in,
					thick);
			}
		}
		if(ans=='6')
		{
			printf("give minimum and maximum latitude\n");
			scanf("%lf%lf",&lat1,&lat2);
			printf("give minimum and maximum longitude\n");
			scanf("%lf%lf",&lon1,&lon2);
			printf("give sample rate\n");
			scanf("%lf",&sample_rate);
			printf("Give output file name.\n");
			scanf("%s",string);
			for(k=0;k<strlen(string);k++)
				if(string[k]=='.')
					string[k]='\0';
			if(scan_height_sample(lat1,lat2,lon1,lon2,&max,&min,sample_rate)>=0)
			{
				printf("Elevation:  min = %ld  max = %ld meters\n",min,max);
				sample(lat1,lat2,lon1,lon2,(double)max,(double)min,
					sample_rate,string);
			}
			else
				printf("scan_height failed.\n\n");
		}
		if(ans=='4')
		{
			printf("give minimum and maximum latitude\n");
			scanf("%lf%lf",&lat1,&lat2);
			printf("give minimum and maximum longitude\n");
			scanf("%lf%lf",&lon1,&lon2);
			printf("Give elevation scale exageration\n");
			scanf("%lf",&r);
			clat=(lat2+lat1)/2.0;
			printf("Give latitude and longitude of sub boxes.\n\n");
			scanf("%lf%lf",&dlat,&dlon);
			row=(lat2-lat1)/dlat+0.5;
			col=(lon2-lon1)/dlon+0.5;
			printf("Working %d rows and %d columns of boxes.\n\n",row,col);
			if(scan_height(lat1,lat2,lon1,lon2,&max,&min)>=0)
			{
				printf("min = %ld  max = %ld\n",min,max);
				km_in=9.0/(dlat*111.7);
				for(i=0;i<row;i++)
				{
					lt2=lat2-dlat*(double)i;
					lt1=lat2-dlat*(double)(i+1);
					for(j=0;j<col;j++)
					{
						ln1=lon1+dlon*(double)j;
						ln2=lon1+dlon*(double)(j+1);
						sprintf(string,"elev%02d%02d.wvf",i+1,j+1);
						printf("%cworking '%s'",cr,string);
						wf_flat(lt1,lt2,ln1,ln2,clat,r,min,max,string,km_in,
							(double)0.125);
					}
				}
			}
		}
		if(ans=='5')	/* projected skin set */
		{
			FpTxt=fopen("elev.txt","wt");
			if(!FpTxt)
			{
				printf("Could not open 'elev.txt' to write prjection data.\n\n");
				printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
				printf("etopo5 l:\\3d\\\n\n");
				exit(0);
			}
			printf("give minimum and maximum latitude (whole degrees)\n");
			scanf("%lf%lf",&lat1,&lat2);
			printf("give minimum and maximum longitude (whole degrees)\n");
			scanf("%lf%lf",&lon1,&lon2);
			printf("Give vertical (elevation) scale exageration\n");
			scanf("%lf",&r);
			clat=(lat2+lat1)/2.0;
		printf("Give latitude and longitude of sub boxes (whole degrees).\n\n");
			scanf("%lf%lf",&dlat,&dlon);
			row=(lat2-lat1)/dlat+0.5;
			col=(lon2-lon1)/dlon+0.5;
			printf("give thickness in inches.\n");
			scanf("%lf",&thick);
			printf("Working %d rows and %d columns of boxes.\n\n",row,col);
			if(choose_projection(&Map,FpTxt)<0)
			{
				printf("No projection chosen.\n\n");
				exit(0);
			}
			get_constants(&Map,FpTxt);
			if(FpTxt)
				fclose(FpTxt);			
			printf("\nscanning for max/min elevation ...\n");
			if(scan_height(lat1,lat2,lon1,lon2,&max,&min)>=0)
			{
				printf("min = %ld  max = %ld\n",min,max);
				maxdely=maxdelx=0.0;
				lon=Map.lon_0/DEG_RAD;
				for(i=0;i<row;i++)
				{
					lt2=lat2-dlat*(double)i;
					lt1=lat2-dlat*(double)(i+1);
					ll2xy(lt1,lon,&xmax,&ymax,Map);
					xmin=xmax;
					ymin=ymax;
					ll2xy(lt1,lon+dlon/2.0,&x,&y,Map);
					maxmin(&xmax,&xmin,x,&ymax,&ymin,y);
					ll2xy(lt1,lon-dlon/2.0,&x,&y,Map);
					maxmin(&xmax,&xmin,x,&ymax,&ymin,y);
					ll2xy(lt2,lon,&x,&y,Map);
					maxmin(&xmax,&xmin,x,&ymax,&ymin,y);
					ll2xy(lt2,lon+dlon/2.0,&x,&y,Map);
					maxmin(&xmax,&xmin,x,&ymax,&ymin,y);
					ll2xy(lt2,lon-dlon/2.0,&x,&y,Map);
					maxmin(&xmax,&xmin,x,&ymax,&ymin,y);
					if(xmax-xmin>maxdelx)
						maxdelx=xmax-xmin;
					if(ymax-ymin>maxdely)
						maxdely=ymax-ymin;
				}
				km_in=maxdely/9.0;
				if(maxdelx>maxdely)
					km_in=maxdelx/9.0;
				for(i=0;i<row;i++)
				{
					lt2=lat2-dlat*(double)i;
					lt1=lat2-dlat*(double)(i+1);
					for(j=0;j<col;j++)
					{
						ln1=lon1+dlon*(double)j;
						ln2=lon1+dlon*(double)(j+1);
						sprintf(string,"elev%02d%02d.wvf",i+1,j+1);
						printf("%cworking '%s'",cr,string);
						wf_flat_skin_proj(lt1,lt2,ln1,ln2,clat,r,
							min,max,string,km_in,thick);
					}
				}
			}
			exit(0);
		}
		if(ans=='a')
		{
			printf("give minimum and maximum latitude\n");
			scanf("%lf%lf",&lat1,&lat2);
			printf("give minimum and maximum longitude\n");
			scanf("%lf%lf",&lon1,&lon2);
			if(scan_height(lat1,lat2,lon1,lon2,&max,&min)>=0)
				printf("min = %ld  max = %ld\n",min,max);
		}
	}while(ans!='x');
}


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

int sample(double lat1,double lat2,double lon1,double lon2,
		double maxz,double minz,double sample_rate,char *name)

{
	int i,j,k;
	int i0,j0,i1,j1,row,col,isamp=sample_rate;
	int arow,acol;
	char string[100],cr=0xd;
	FILE *fplbs,*fpin,*fpout;
	double dellat,dellon,sign180;
	long drow=ROW,dcol=COL,offset;


	dellon=dellat=sample_rate/(double)12.0;
	if(lat1==-90.0)
		lat1=-90.0+1.0/12.0;
	i0=(90.0-lat2)*12.0;
	i1=(90.0-lat1)*12.0;
	if(lon1<0.0)
	{
		lon1+=360.0;
		sign180=-1;
	}
	j0=lon1*(double)12.0;
	if(lon2<0.0)
	{
		lon2+=360.0;
		sign180=-1;
	}
	if(lon2==0.0||lon2==360.0)
		lon2=360.0-1.0/12.0;
	if(lat1>lat2)
	{
		printf("Minimum latitude is larger than maximum.\n\n");
		exit(0);
	}
	if(lon1>lon2)
	{
		printf("Minimum longitude is larger than maximum.\n\n");
		exit(0);
	}
	j1=lon2*(double)12.0;
	row=i1-i0+1;
	col=j1-j0+1;
	arow=(row+isamp-1)/isamp;
	acol=(col+isamp-1)/isamp;
	lat1=lat2-(double)(arow-1)*dellat;
	lon2=lon1+(double)(acol-1)*dellon;
	lat2=(double)90.0-(double)i0/12.0;
	lat1=(double)90.0-(double)i1/12.0;
	lon1=(double)j0/(double)12.0;
	lon2=(double)j1/(double)12.0;
	if(lon1>=180.0)
		lon1-=360.0;
	if(lon2>180.0)
		lon2-=360.0;
	if(lon1>=0.0&&lon2<0.0||lon1>lon2)
		lon2+=360.0;
	sprintf(string,"%s.lbs",name);
	fplbs=fopen(string,"wt");
	if(!fplbs)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	else
	{
		fprintf(fplbs,"FILE_TYPE          = ELEVATION_ARRAY\n");
		fprintf(fplbs,"IMAGE_LINES        = %6d\n",arow);
		fprintf(fplbs,"LINE_SAMPLES       = %6d\n",acol);
		fprintf(fplbs,"UPPER_LEFT_LAT     = %10.5lf\n",lat2);
		fprintf(fplbs,"UPPER_LEFT_LON     = %10.5lf\n",lon1);
		fprintf(fplbs,"LOWER_RIGHT_LAT    = %10.5lf\n",lat1);
		fprintf(fplbs,"LOWER_RIGHT_LON    = %10.5lf\n",lon2);
		fprintf(fplbs,"DELTA_LAT          = %10.6lf\n",dellat);
		fprintf(fplbs,"DELTA_LON          = %10.6lf\n",dellon);
		fprintf(fplbs,"NUMBER_OF_LAT      = %6d\n",arow);
		fprintf(fplbs,"NUMBER_OF_LON      = %6d\n",acol);
		fprintf(fplbs,"MAXIMUM_Z          = %15.5lf\n",maxz);
		fprintf(fplbs,"MINIMUM_Z          = %15.5lf\n",minz);
		fprintf(fplbs,"SOURCE_OF_DATA     = ETOPO5\n");
		fprintf(fplbs,"DATA_TYPE          = short\n");
		fprintf(fplbs,"DATA_SIZE          = 2 bytes\n");
		fprintf(fplbs,"DATA_FILE_POINTER  = '%s.D16'\n",name);
	}
	fclose(fplbs);
	sprintf(string,"%s.d16",name);
	fpout=fopen(string,"wb");
	if(!fpout)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	offset=(i0*dcol+(long)j0)*2;
	fpin=fopen_path("globe16.dat","rb");
	if(!fpin)
	{
		printf("Could not open 'globe16.dat' to read.\n\n");
		printf("If file is not on this disk include path, e.g.\n");
		printf("\netopo5 l:\\3d\\\n");
		exit(0);
	}
	fseek(fpin,offset,SEEK_SET);
	for(i=0;i<row;i++)
	{
		if(i%isamp==0)
		{
			fread((char *)Elev,sizeof(int),col,fpin);
			printf("%c%d",cr,i/isamp+1);
			for(j=0;j<acol;j++)
				E16[j]=Elev[j*isamp];
			fwrite((char *)E16,sizeof(int),acol,fpout);
		}
		offset+=dcol*2;
		fseek(fpin,offset,SEEK_SET);
	}
	printf("\n");
	fclose(fpin);
	fclose(fpout);
}


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

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);
}


/**********************************************************************
**
**	thickness is in inches
**
**********************************************************************/

int wf_flat_skin_proj(double lat1,double lat2,double lon1,double lon2,
		double clat,double r,long dmin,long dmax,char *name,double km_in,
		double thick)

{
	int i,j,k,l,m,n;
	long drow=ROW,dcol=COL;
	int i0,j0,i1,j1,row,col;
	double lat,lon,x,y,z,rho,ln1=lon1,ln2=lon2;
	long offset,pts;
	long ij[2][2][2];	/* cube [depth][col][row] */
	FILE *fpin,*fpout,*fplbl,*fpdat,*fptxt;
	long drange,dval,dnum=0;
	char string[100],cr=0xd;
	double scale=cos(clat*DEG_RAD);
	long nm,cl,rw,lr,lc,lrow,lcol,l2=2,l1=1,l3=3,l4=4;
	double mlon,mlat,mz;
	double mx,my;
	int sign180=1;


	drange=dmax-dmin;
	i0=(90.0-lat2)*12.0;
	i1=(90.0-lat1)*12.0;
	if(lon1<0.0)
	{
		lon1+=360.0;
		sign180=-1;
	}
	j0=lon1*(double)12.0;
	if(lon2<0.0)
	{
		lon2+=360.0;
		sign180=-1;
	}
	j1=lon2*(double)12.0;
	row=i1-i0+1;
	col=j1-j0+1;
	sprintf(string,"elev.lbl");
	fplbl=fopen(string,"wt");
	if(!fplbl)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	else
	{
		fprintf(fplbl,"FILE_TYPE          = IMAGE\n");
		fprintf(fplbl,"IMAGE_LINES        =  %4d\n",row);
		fprintf(fplbl,"LINE_SAMPLES       =  %4d\n",col);
		fprintf(fplbl,"IMAGE_POINTER      = 'elev.dat'\n");
		fclose(fplbl);
	}
	fclose(fplbl);
	sprintf(string,"elev.dat");
	fpdat=fopen(string,"wb");
	if(!fpdat)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	offset=(i0*dcol+(long)j0)*2;
	lat2=(double)90.0-(double)i0/12.0;
	lat1=(double)90.0-(double)i1/12.0;
	lon1=(double)j0/(double)12.0;
	lon2=(double)j1/(double)12.0;
	mlon=(lon1+lon2)/2.0;
	if(mlon>180.0)
		mlon-=360;
	mlat=(lat1+lat2)/2.0;
	ll2xy(mlat,mlon,&mx,&my,Map);
	mz=(double)(dmax+dmin)/2.0/1000.0*r/km_in;
	fpin=fopen_path("globe16.dat","rb");
	if(!fpin)
	{
		printf("Could not open 'globe16.dat' to read.\n\n");
		printf("If file is not on this disk include path, e.g.\n");
		printf("\netopo5 l:\\3d\\\n");
		exit(0);
	}
	fpout=fopen(name,"wt");
	if(!fpout)
	{
		printf("Could not open '%s' to write.\n\n",name);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	fprintf(fpout,"# lat = %lf to %lf\n",lat1,lat2);
	fprintf(fpout,"# lon = %lf to %lf\n",ln1,ln2);
	fprintf(fpout,"# central lat = %lf\n",clat);
	fprintf(fpout,"# vertical exageration = %lf\n",r);
	fprintf(fpout,"# maximum elevation = %7ld.0 m\n",dmax);
	fprintf(fpout,"#  middle elevation = %9.1lf m\n",mz*km_in/r*1000.0);
	fprintf(fpout,"# minimum elevation = %7ld.0 m\n",dmin);
	fprintf(fpout,"# horizontal scale = %7.3lf km/in\n",km_in);
	fprintf(fpout,"# flat skin model -- projected\n");
	fptxt=fopen_path("elev.txt","rt");
	if(fptxt)
	{
		while(fgets(string,100,fptxt))
		{
			fprintf(fpout,"# %s",string);
		}
		fclose(fptxt);
	}
	fprintf(fpout,"# skin thickness = %6.3lf in\n",thick);
	fprintf(fpout,"# array size = %d rows x %d cols\n",row,col);
	fprintf(fpout,"# %ld verticies\n",(long)row*(long)col*(long)2);
	fprintf(fpout,"# %ld triangles\n",(long)4*((long)row*(long)col-(long)1));
	fprintf(fpout,"g skin\n");
	fseek(fpin,offset,SEEK_SET);
	for(i=0;i<row;i++)
	{
		printf("%c%d",cr,i+1);
		lat=(lat2-(double)i/12.0);
		fread((char *)Elev,sizeof(int),col,fpin);
		offset+=dcol*2;
		fseek(fpin,offset,SEEK_SET);
		for(j=0;j<col;j++)
		{
			lon=(lon1+(double)j/12.0);
			if(sign180>0&&lon>180.0)
				lon-=360.0;
			if(sign180<0&&lon>=180.0)
				lon-=360.0;
			ll2xy(lat,lon,&x,&y,Map);
			rho=Elev[j];
			rho=rho/1000.0;
			x=(x-mx)/km_in;
			y=(y-my)/km_in;
			z=rho*r/km_in-mz;
			fprintf(fpout,"v %12.6lf %12.6lf %12.6lf\n",x,y,z);
			fprintf(fpout,"v %12.6lf %12.6lf %12.6lf\n",x,y,z-thick);
			dval=Elev[j];
			dval-=dmin;
			dval*=255;
			E8[j]=dval/drange;
		}
		fwrite((char *)E8,sizeof(char),col,fpdat);
	}
	printf("\n");
	fprintf(fpout,"# %ld triangles\n",(long)4*((long)row*(long)col-(long)1));
	fprintf(fpout,"# lat = %lf to %lf\n",lat1,lat2);
	fprintf(fpout,"# lon = %lf to %lf\n",lon1,lon2);
	fprintf(fpout,"# flat skin model -- projected\n");
	fprintf(fpout,"g skin\n");
	lrow=row;
	lcol=col;
	for(i=1;i<row;i++)
	{
		lr=i;
		ij[0][0][0]=lcol*l2*(lr-l1)+l1;	/* top upper left (1st cube) */
		ij[0][1][0]=lcol*l2*(lr-l1)+l3;	/* top upper right (1st cube) */
		ij[0][0][1]=lcol*l2*lr+l1;	/* top lower left (1st cube) */
		ij[0][1][1]=lcol*l2*lr+l3;	/* top lower right (1st cube) */

		ij[1][0][0]=lcol*l2*(lr-l1)+l2;	/* bottom upper left (1st cube) */
		ij[1][1][0]=lcol*l2*(lr-l1)+l4;	/* bottom upper right (1st cube) */
		ij[1][0][1]=lcol*l2*lr+l2;	/* bottom lower left (1st cube) */
		ij[1][1][1]=lcol*l2*lr+l4;	/* bottom lower right (1st cube) */
		for(j=1;j<col;j++)
		{
			lc=j;
/* top surface */
			fprintf(fpout,"f %ld %ld %ld\n",ij[0][0][1],ij[0][1][0],ij[0][0][0]);
			fprintf(fpout,"f %ld %ld %ld\n",ij[0][0][1],ij[0][1][1],ij[0][1][0]);
/* bottom surface */
			fprintf(fpout,"f %ld %ld %ld\n",ij[1][0][1],ij[1][0][0],ij[1][1][0]);
			fprintf(fpout,"f %ld %ld %ld\n",ij[1][0][1],ij[1][1][0],ij[1][1][1]);
/* left edge */
			if(j==1)
			{
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][0][0],ij[0][0][1],ij[0][0][0]);
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][0][0],ij[1][0][1],ij[0][0][1]);
			}
/* right edge */
			if(j==col-1)
			{
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][1][1],ij[0][1][0],ij[0][1][1]);
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][1][1],ij[1][1][0],ij[0][1][0]);
			}
/* upper edge */
			if(i==1)
			{
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][1][0],ij[1][0][0],ij[0][0][0]);
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][1][0],ij[0][0][0],ij[0][1][0]);
			}
/* lower edge */
			if(i==row-1)
			{
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][0][1],ij[0][1][1],ij[0][0][1]);
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][0][1],ij[1][1][1],ij[0][1][1]);
			}
			for(l=0;l<2;l++)	/* move cube 1 to left */
				for(m=0;m<2;m++)
					for(n=0;n<2;n++)
						ij[l][m][n]+=2;
		}
	}
	j=col;
	j=0;
	fclose(fpin);
	fclose(fpout);
	fclose(fpdat);
}




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

int maxmin(double *xmax,double *xmin,double x,
			  double *ymax,double *ymin,double y)

{
	if(x>*xmax)
		*xmax=x;
	if(x<*xmin)
		*xmin=x;
	if(y>*ymax)
		*ymax=y;
	if(y<*ymin)
		*ymin=y;
}






/**********************************************************************
**
**	thickness is in inches
**
**********************************************************************/

int wf_flat_skin(double lat1,double lat2,double lon1,double lon2,
		double clat,double r,long dmin,long dmax,char *name,double km_in,
		double thick)

{
	int i,j,k,l,m,n;
	long drow=ROW,dcol=COL;
	int i0,j0,i1,j1,row,col;
	double lat,lon,x,y,z,rho,ln1=lon1,ln2=lon2;
	long offset,pts;
	long ij[2][2][2];	/* cube [depth][col][row] */
	FILE *fpin,*fpout,*fplbl,*fpdat;
	long drange,dval,dnum=0;
	char string[100],cr=0xd;
	double scale=cos(clat*DEG_RAD);
	long nm,cl,rw,lr,lc,lrow,lcol,l2=2,l1=1,l3=3,l4=4;
	double mlon,mlat,mz;


	drange=dmax-dmin;
	i0=(90.0-lat2)*12.0;
	i1=(90.0-lat1)*12.0;
	if(lon1<0.0)
		lon1+=360.0;
	j0=lon1*(double)12.0;
	if(lon2<0.0)
		lon2+=360.0;
	j1=lon2*(double)12.0;
	row=i1-i0+1;
	col=j1-j0+1;
	sprintf(string,"elev.lbl");
	fplbl=fopen(string,"wt");
	if(!fplbl)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	else
	{
		fprintf(fplbl,"FILE_TYPE          = IMAGE\n");
		fprintf(fplbl,"IMAGE_LINES        =  %4d\n",row);
		fprintf(fplbl,"LINE_SAMPLES       =  %4d\n",col);
		fprintf(fplbl,"IMAGE_POINTER      = 'elev.dat'\n");
		fclose(fplbl);
	}
	fclose(fplbl);
	sprintf(string,"elev.dat");
	fpdat=fopen(string,"wb");
	if(!fpdat)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}

	offset=(i0*dcol+(long)j0)*2;
	lat2=(double)90.0-(double)i0/12.0;
	lat1=(double)90.0-(double)i1/12.0;
	lon1=(double)j0/(double)12.0;
	lon2=(double)j1/(double)12.0;
	mlon=(lon1+lon2)/2.0;
	mlat=(lat1+lat2)/2.0;
	mz=(double)(dmax+dmin)/2.0/1000.0*r*km_in;
	fpin=fopen_path("globe16.dat","rb");
	if(!fpin)
	{
		printf("Could not open 'globe16.dat' to read.\n\n");
		printf("If file is not on this disk include path, e.g.\n");
		printf("\netopo5 l:\\3d\\\n");
		exit(0);
	}
	fpout=fopen(name,"wt");
	if(!fpout)
	{
		printf("Could not open '%s' to write.\n\n",name);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	fprintf(fpout,"# lat = %lf to %lf\n",lat1,lat2);
	fprintf(fpout,"# lon = %lf to %lf\n",ln1,ln2);
	fprintf(fpout,"# central lat = %lf\n",clat);
	fprintf(fpout,"# vertical exageration = %lf\n",r);
	fprintf(fpout,"# horizontal scale = %7.3lf km/in\n",(double)1.0/km_in);
	fprintf(fpout,"# flat skin model\n");
	fprintf(fpout,"# skin thickness = %6.3lf in\n",thick);
	fprintf(fpout,"# array size = %d rows x %d cols\n",row,col);
	fprintf(fpout,"# %ld verticies\n",(long)row*(long)col*(long)2);
	fprintf(fpout,"# %ld triangles\n",(long)4*((long)row*(long)col-(long)1));
	fprintf(fpout,"g skin\n");
	fseek(fpin,offset,SEEK_SET);
	for(i=0;i<row;i++)
	{
		printf("%c%d",cr,i+1);
		lat=(lat2-(double)i/12.0);
		y=-(mlat-lat)/(double)360.0*E_CIR*km_in;
		fread((char *)Elev,sizeof(int),col,fpin);
		offset+=dcol*2;
		fseek(fpin,offset,SEEK_SET);
		for(j=0;j<col;j++)
		{
			lon=(lon1+(double)j/12.0);
			rho=Elev[j];
			rho=rho/1000.0;
			x=(lon-mlon)/(double)360.0*E_CIR*scale*km_in;
			z=rho*r*km_in-mz;
			fprintf(fpout,"v %12.6lf %12.6lf %12.6lf\n",x,y,z);
			fprintf(fpout,"v %12.6lf %12.6lf %12.6lf\n",x,y,z-thick);
			dval=Elev[j];
			dval-=dmin;
			dval*=255;
			E8[j]=dval/drange;
		}
		fwrite((char *)E8,sizeof(char),col,fpdat);
	}
	printf("\n");
	fprintf(fpout,"# %ld triangles\n",(long)4*((long)row*(long)col-(long)1));
	fprintf(fpout,"# lat = %lf to %lf\n",lat1,lat2);
	fprintf(fpout,"# lon = %lf to %lf\n",lon1,lon2);
	fprintf(fpout,"# flat skin model\n");
	fprintf(fpout,"g skin\n");
	lrow=row;
	lcol=col;
	for(i=1;i<row;i++)
	{
		lr=i;
		ij[0][0][0]=lcol*l2*(lr-l1)+l1;	/* top upper left (1st cube) */
		ij[0][1][0]=lcol*l2*(lr-l1)+l3;	/* top upper right (1st cube) */
		ij[0][0][1]=lcol*l2*lr+l1;	/* top lower left (1st cube) */
		ij[0][1][1]=lcol*l2*lr+l3;	/* top lower right (1st cube) */

		ij[1][0][0]=lcol*l2*(lr-l1)+l2;	/* bottom upper left (1st cube) */
		ij[1][1][0]=lcol*l2*(lr-l1)+l4;	/* bottom upper right (1st cube) */
		ij[1][0][1]=lcol*l2*lr+l2;	/* bottom lower left (1st cube) */
		ij[1][1][1]=lcol*l2*lr+l4;	/* bottom lower right (1st cube) */
		for(j=1;j<col;j++)
		{
			lc=j;
/* top surface */
			fprintf(fpout,"f %ld %ld %ld\n",ij[0][0][1],ij[0][1][0],ij[0][0][0]);
			fprintf(fpout,"f %ld %ld %ld\n",ij[0][0][1],ij[0][1][1],ij[0][1][0]);
/* bottom surface */
			fprintf(fpout,"f %ld %ld %ld\n",ij[1][0][1],ij[1][0][0],ij[1][1][0]);
			fprintf(fpout,"f %ld %ld %ld\n",ij[1][0][1],ij[1][1][0],ij[1][1][1]);
/* left edge */
			if(j==1)
			{
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][0][0],ij[0][0][1],ij[0][0][0]);
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][0][0],ij[1][0][1],ij[0][0][1]);
			}
/* right edge */
			if(j==col-1)
			{
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][1][1],ij[0][1][0],ij[0][1][1]);
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][1][1],ij[1][1][0],ij[0][1][0]);
			}
/* upper edge */
			if(i==1)
			{
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][1][0],ij[1][0][0],ij[0][0][0]);
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][1][0],ij[0][0][0],ij[0][1][0]);
			}
/* lower edge */
			if(i==row-1)
			{
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][0][1],ij[0][1][1],ij[0][0][1]);
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][0][1],ij[1][1][1],ij[0][1][1]);
			}
			for(l=0;l<2;l++)	/* move cube 1 to left */
				for(m=0;m<2;m++)
					for(n=0;n<2;n++)
						ij[l][m][n]+=2;
		}
	}
	j=col;
	j=0;
	fclose(fpin);
	fclose(fpout);
	fclose(fpdat);
}



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

int scan_height(double lat1,double lat2,double lon1,double lon2,
		long *max,long *min)

{
	int i,j,k,m,n;
	long drow=ROW,dcol=COL;
	int i0,j0,i1,j1,row,col;
	double lat,lon,x,y,z,rho;
	long offset,pts,ij[2][2];
	FILE *fpin;
	char string[100],cr=0xd;
	long nm,cl,rw,le;

	if(lon2==0.0||lon2==360.0)
		lon2=360.0-1.0/12.0;
	if(lat1==-90.0)
		lat1=-90.0+1.0/12.0;
	*max=-1000000000,*min=1000000000;
	i0=(90.0-lat2)*12.0;
	if(i0<0||i0>=ROW)
	{
		printf("maximum latitude is not valid.\n\n");
		printf("lat_2 = %lf -- i0 = %d must be between 0 and 2159.\n\n",
			lat2,i0);
		return(-1);
	}
	i1=(90.0-lat1)*12.0;
	if(i1<0||i1>=ROW)
	{
		printf("minimum latitude is not valid.\n\n");
		printf("lat_1 = %lf -- i1 = %d must be between 0 and 2159.\n\n",
			lat1,i1);
		return(-1);
	}
	if(lon1<0.0)
		lon1+=360.0;
	j0=lon1*(double)12.0;
	if(j0<0||j0>=COL)
	{
		printf("minimum longitude is not valid.\n\n");
		printf("lon_1 = %lf -- j0 = %d must be between 0 and 4319.\n\n",
			lon1,j0);
		return(-1);
	}
	if(lon2<0.0)
		lon2+=360.0;
	j1=lon2*(double)12.0;
	if(j1<0||j1>=COL)
	{
		printf("maximum longitude is not valid.\n\n");
		printf("lon_2 = %lf -- j1 = %d must be between 0 and 4319.\n\n",
			lon2,j1);
		return(-1);
	}
	row=i1-i0+1;
	col=j1-j0+1;
	if(row<0||row>ROW)
	{
		printf("Latitude range is not valid.\n\n");
		printf("row = %d - %d + 1 = %d\n",i1,i0,row);
		return(-1);
	}
	if(col<0||col>COL)
	{
		printf("longitude range is not valid.\n\n");
		printf("col = %d - %d + 1 = %d\n",j1,j0,col);
		printf(
			"You may have crossed the 0 deg meridian which is not allowed.\n\n");
		return(-1);
	}
	printf("Processing %d rows and %d columns\n",row,col);
	pts=(long)row*(long)col+(long)((row-1)*2)+(long)((col-1)*2);
	offset=(i0*dcol+(long)j0)*2;
	lat2=(double)90.0-(double)i0/12.0;
	lat1=(double)90.0-(double)i1/12.0;
	lon1=(double)j0/(double)12.0;
	lon2=(double)j1/(double)12.0;
	fpin=fopen_path("globe16.dat","rb");
	if(!fpin)
	{
		printf("Could not open 'globe16.dat' to read.\n\n");
		printf("If file is not on this disk include path, e.g.\n");
		printf("\netopo5 l:\\3d\\\n");
		exit(0);
	}
	fseek(fpin,offset,SEEK_SET);
	for(i=0;i<row;i++)
	{
		printf("%c%5d of %d",cr,i+1,row);
		lat=(lat2-(double)i/12.0);
		y=-(lat2-lat)/(double)360.0*E_CIR;
		fread((char *)Elev,sizeof(int),col,fpin);
		offset+=dcol*2;
		fseek(fpin,offset,SEEK_SET);
		for(j=0;j<col;j++)
		{
			le=Elev[j];
			if(le>*max)
				*max=le;
			if(le<*min)
				*min=le;
		}
	}
	printf("\n");
	if(fpin)
		fclose(fpin);
	return(1);
}

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

int scan_height_sample(double lat1,double lat2,double lon1,double lon2,
		long *max,long *min,double sample)

{
	int i,j,k,m,n;
	long drow=ROW,dcol=COL;
	int i0,j0,i1,j1,row,col,isamp;
	double lat,lon,x,y,z,rho;
	long offset,pts,ij[2][2];
	FILE *fpin;
	char string[100],cr=0xd;
	long nm,cl,rw,le;

	isamp=sample;
	*max=-1000000000,*min=1000000000;
	if(lat1==-90.0)
		lat1=-90.0+1.0/12.0;
	i0=(90.0-lat2)*12.0;
	if(lon2==0.0||lon2==360.0)
		lon2=360.0-1.0/12.0;
	if(i0<0||i0>=ROW)
	{
		printf("maximum latitude is not valid.\n\n");
		printf("lat_2 = %lf -- i0 = %d must be between 0 and 2159.\n\n",
			lat2,i0);
		return(-1);
	}
	i1=(90.0-lat1)*12.0;
	if(i1<0||i1>=ROW)
	{
		printf("minimum latitude is not valid.\n\n");
		printf("lat_1 = %lf -- i1 = %d must be between 0 and 2159.\n\n",
			lat1,i1);
		return(-1);
	}
	if(lon1<0.0)
		lon1+=360.0;
	j0=lon1*(double)12.0;
	if(j0<0||j0>=COL)
	{
		printf("minimum longitude is not valid.\n\n");
		printf("lon_1 = %lf -- j0 = %d must be between 0 and 4319.\n\n",
			lon1,j0);
		return(-1);
	}
	if(lon2<0.0)
		lon2+=360.0;
	j1=lon2*(double)12.0;
	if(j1<0||j1>=COL)
	{
		printf("maximum longitude is not valid.\n\n");
		printf("lon_2 = %lf -- j1 = %d must be between 0 and 4319.\n\n",
			lon2,j1);
		return(-1);
	}
	row=i1-i0+1;
	col=j1-j0+1;
	if(row<0||row>ROW)
	{
		printf("Latitude range is not valid.\n\n");
		printf("row = %d - %d + 1 = %d\n",i1,i0,row);
		return(-1);
	}
	if(col<0||col>COL)
	{
		printf("longitude range is not valid.\n\n");
		printf("col = %d - %d + 1 = %d\n",j1,j0,col);
		printf(
			"You may have crossed the 0 deg meridian which is not allowed.\n\n");
		return(-1);
	}
	printf("Processing %d rows and %d columns to find max/min elevation.\n",
		row/isamp+1,col/isamp+1);
	pts=(long)row*(long)col+(long)((row-1)*2)+(long)((col-1)*2);
	offset=(i0*dcol+(long)j0)*2;
	lat2=(double)90.0-(double)i0/12.0;
	lat1=(double)90.0-(double)i1/12.0;
	lon1=(double)j0/(double)12.0;
	lon2=(double)j1/(double)12.0;
	fpin=fopen_path("globe16.dat","rb");
	if(!fpin)
	{
		printf("Could not open 'globe16.dat' to read.\n\n");
		printf("If file is not on this disk include path, e.g.\n");
		printf("\netopo5 l:\\3d\\\n");
		exit(0);
	}
	fseek(fpin,offset,SEEK_SET);
	for(i=0;i<row;i++)
	{
		lat=(lat2-(double)i/12.0);
		y=-(lat2-lat)/(double)360.0*E_CIR;
		if(i%isamp==0)
		{
			fread((char *)Elev,sizeof(int),col,fpin);
			printf("%c%5d of %d rows",cr,i/isamp+1,row/isamp+1);
			for(j=0;j<col;j++)
			{
				if(j%isamp==0)
				{
					le=Elev[j];
					if(le>*max)
						*max=le;
					if(le<*min)
						*min=le;
				}
			}
		}
		offset+=dcol*2;
		fseek(fpin,offset,SEEK_SET);
	}
	printf("\n");
	if(fpin)
		fclose(fpin);
	return(1);
}




/**********************************************************************
**
**	bottom thickness is in inches
**
**********************************************************************/

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 i,j,k,m,n;
	long drow=ROW,dcol=COL;
	int i0,j0,i1,j1,row,col;
	double lat,lon,x,y,z,rho,ln1=lon1,ln2=lon2,mlat,mlon;
	long offset,pts,ij[2][2],dmid;
	FILE *fpin,*fpout,*fplbl,*fpdat;
	long drange,dval,dnum=0;
	char string[100];
	double scale=cos(clat*DEG_RAD);
	long nm,cl,rw;

	drange=dmax-dmin;
	i0=(90.0-lat2)*12.0;
	i1=(90.0-lat1)*12.0;
	if(lon1<0.0)
		lon1+=360.0;
	j0=lon1*(double)12.0;
	if(lon2<0.0)
		lon2+=360.0;
	j1=lon2*(double)12.0;
	mlat=(lat1+lat2)/2.0;
	mlon=(lon1+lon2)/2.0;
	dmid=(dmin+dmax)/2.0;
	row=i1-i0+1;
	col=j1-j0+1;
	sprintf(string,"elev.lbl");
	fplbl=fopen(string,"wt");
	if(!fplbl)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	else
	{
		fprintf(fplbl,"FILE_TYPE          = IMAGE\n");
		fprintf(fplbl,"IMAGE_LINES        =  %4d\n",row);
		fprintf(fplbl,"LINE_SAMPLES       =  %4d\n",col);
		fprintf(fplbl,"IMAGE_POINTER      = 'elev.dat'\n");
		fclose(fplbl);
	}
	fclose(fplbl);
	sprintf(string,"elev.dat");
	fpdat=fopen(string,"wb");
	if(!fpdat)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}

	pts=(long)row*(long)col+(long)((row-1)*2)+(long)((col-1)*2);
	offset=(i0*dcol+(long)j0)*2;
	lat2=(double)90.0-(double)i0/12.0;
	lat1=(double)90.0-(double)i1/12.0;
	lon1=(double)j0/(double)12.0;
	lon2=(double)j1/(double)12.0;
	fpin=fopen_path("globe16.dat","rb");
	if(!fpin)
	{
		printf("Could not open 'globe16.dat' to read.\n\n");
		printf("If file is not on this disk include path, e.g.\n");
		printf("\netopo5 l:\\3d\\\n");
		exit(0);
	}
	fpout=fopen(name,"wt");
	if(!fpout)
	{
		printf("Could not open '%s' to write.\n\n",name);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	fprintf(fpout,"# lat = %lf to %lf\n",lat1,lat2);
	fprintf(fpout,"# lon = %lf to %lf\n",ln1,ln2);
	fprintf(fpout,"# central lat = %lf\n",clat);
	fprintf(fpout,"# mid lat = %lf\n",mlat);
	fprintf(fpout,"# min lon = %lf\n",mlon);
	fprintf(fpout,"# mid z   = %l\n",dmid);
	fprintf(fpout,"# vertical exageration = %lf\n",r);
	fprintf(fpout,"# %7.3lf km/in\n",(double)1.0/km_in);
	fprintf(fpout,"# %6.3lf in bottom thickness.\n",bottom_thick);
	fprintf(fpout,"g globe\n");
	fprintf(fpout,"# %d rows x %d cols\n",row,col);
	fprintf(fpout,"g\n");
	fseek(fpin,offset,SEEK_SET);
	for(i=0;i<row;i++)
	{
		lat=(lat2-(double)i/12.0);
		y=-(mlat-lat)/(double)360.0*E_CIR*km_in;
		fread((char *)Elev,sizeof(int),col,fpin);
		offset+=dcol*2;
		fseek(fpin,offset,SEEK_SET);
		for(j=0;j<col;j++)
		{
			lon=(lon1+(double)j/12.0);
			rho=(long)Elev[j]-dmid;
			rho=rho/1000.0;
			x=(lon-mlon)/(double)360.0*E_CIR*scale*km_in;
			z=rho*r*km_in;
			fprintf(fpout,"v %12.6lf %12.6lf %12.6lf\n",x,y,z);
			dval=Elev[j];
			dval-=dmin;
			dval*=255;
			E8[j]=dval/drange;
		}
		fwrite((char *)E8,sizeof(char),col,fpdat);
	}
	z=dmin;
	z=z/1000.0*r*km_in-bottom_thick;
	nm=(long)row*(long)col+1;
	cl=col-1;
	rw=row-1;
	i=0;
	lat=(lat2-(double)i/12.0);
	y=-(lat2-lat)/(double)360.0*E_CIR*km_in;
	for(j=0;j<cl;j++)
	{
		lon=(lon1+(double)j/12.0);
		x=(lon-lon1)/(double)360.0*E_CIR*scale*km_in;
		fprintf(fpout,"v %12.6lf %12.6lf %12.6lf\n",x,y,z);
	}
	j=cl;
	lon=(lon1+(double)j/12.0);
	x=(lon-lon1)/(double)360.0*E_CIR*scale*km_in;
	for(i=0;i<rw;i++)
	{
		lat=(lat2-(double)i/12.0);
		y=-(lat2-lat)/(double)360.0*E_CIR*km_in;
		fprintf(fpout,"v %12.6lf %12.6lf %12.6lf\n",x,y,z);
	}
	i=rw;
	lat=(lat2-(double)i/12.0);
	y=-(lat2-lat)/(double)360.0*E_CIR*km_in;
	for(j=cl;j>0;j--)
	{
		lon=(lon1+(double)j/12.0);
		x=(lon-lon1)/(double)360.0*E_CIR*scale*km_in;
		fprintf(fpout,"v %12.6lf %12.6lf %12.6lf\n",x,y,z);
	}
	j=0;
	lon=(lon1+(double)j/12.0);
	x=(lon-lon1)/(double)360.0*E_CIR*scale*km_in;
	for(i=rw;i>0;i--)
	{
		lat=(lat2-(double)i/12.0);
		y=-(lat2-lat)/(double)360.0*E_CIR*km_in;
		fprintf(fpout,"v %12.6lf %12.6lf %12.6lf\n",x,y,z);
	}
	fprintf(fpout,"# %ld vertices\n",pts);
	fprintf(fpout,"# lat = %lf to %lf\n",lat1,lat2);
	fprintf(fpout,"# lon = %lf to %lf\n",lon1,lon2);
	fprintf(fpout,"# earth radius = %lf\n",r);
	fprintf(fpout,"g globe\n");
	for(i=1;i<row;i++)
	{
		ij[1][1]=(long)i*(long)col+2;
		ij[0][0]=(long)(i-1)*(long)col+1;
		ij[0][1]=ij[0][0]+1;
		ij[1][0]=ij[1][1]-1;
		for(j=1;j<col;j++)
		{
			fprintf(fpout,"f %ld %ld %ld\n",ij[1][0],ij[0][1],ij[0][0]);
			fprintf(fpout,"f %ld %ld %ld\n",ij[1][0],ij[1][1],ij[0][1]);
			for(m=0;m<2;m++)
				for(n=0;n<2;n++)
					ij[m][n]+=1;
		}
	}
	i=0;
	for(j=0;j<cl;j++)	/* top -- left to right */
	{
		ij[1][0]=j+1;
		ij[1][1]=j+2;
		ij[0][0]=nm+j;
		ij[0][1]=nm+j+1;
		fprintf(fpout,"f %ld %ld %ld\n",ij[1][0],ij[0][1],ij[0][0]);
		fprintf(fpout,"f %ld %ld %ld\n",ij[1][0],ij[1][1],ij[0][1]);
	}
	j=col;
	for(i=0;i<rw;i++)	/* right -- top to bottom */
	{
		ij[0][0]=(long)(i+1)*(long)col;
		ij[1][0]=(long)(i+2)*(long)col;
		ij[0][1]=nm+cl+i;
		ij[1][1]=nm+cl+i+1;
		fprintf(fpout,"f %ld %ld %ld\n",ij[1][0],ij[0][1],ij[0][0]);
		fprintf(fpout,"f %ld %ld %ld\n",ij[1][0],ij[1][1],ij[0][1]);
	}
	i=row;
	for(j=0;j<cl;j++)	/* bottom -- right to left */
	{
		ij[0][0]=(long)i*(long)col-1-j;
		ij[0][1]=(long)i*(long)col-j;
		ij[1][0]=nm+cl+rw+j+1;
		ij[1][1]=nm+cl+rw+j;
		fprintf(fpout,"f %ld %ld %ld\n",ij[1][0],ij[0][1],ij[0][0]);
		fprintf(fpout,"f %ld %ld %ld\n",ij[1][0],ij[1][1],ij[0][1]);
	}
	j=0;
	for(i=0;i<rw;i++)	/* left -- bottom to top */
	{
		ij[0][1]=(long)(row-i-2)*(long)col+1;
		ij[1][1]=(long)(row-i-1)*(long)col+1;
		ij[0][0]=nm+cl*2+rw+i+1;
		if(ij[0][0]>pts)
			ij[0][0]=nm;
		ij[1][0]=nm+cl*2+rw+i;
		fprintf(fpout,"f %ld %ld %ld\n",ij[1][0],ij[0][1],ij[0][0]);
		fprintf(fpout,"f %ld %ld %ld\n",ij[1][0],ij[1][1],ij[0][1]);
	}
	/* make bottom */
	ij[0][0]=nm;
	ij[0][1]=nm+cl;
	ij[1][1]=nm+cl+rw;
	ij[1][0]=nm+cl*2+rw;
	fprintf(fpout,"f %ld %ld %ld\n",ij[1][0],ij[0][0],ij[0][1]);
	fprintf(fpout,"f %ld %ld %ld\n",ij[1][0],ij[0][1],ij[1][1]);
	pts=(long)(row-1)*(long)(col-1)*2+rw*4+cl*4+2;
	fprintf(fpout,"# %ld elements\n",pts);
	fclose(fpin);
	fclose(fpout);
	fclose(fpdat);
}




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

int	wf(double lat1,double lat2,double lon1,double lon2,double r,
			double zscale,int samp,double thick,char *name)

{
	int i,j,k,m,n;
	long drow=ROW,dcol=COL;
	int i0,j0,i1,j1,row,col;
	long lrow,lcol,mxpts;
	double lat,lon,lon0,x,y,z,rho;
	double xbar=0.0,ybar=0.0,zbar=0.0,scale;
	long offset,pts,ij[2][2];
	long ij3[2][2][2];	/* cube [depth][col][row] */
	long l,lr,lc,l2=2,l1=1,l3=3,l4=4,ppts=0;
	FILE *fpin,*fpout;
	long dmax=9000,dmin=-12000,drange,dval,dnum=0;
	char string[100],cr=0xd;
	int ds=1;

	if(thick<=0.0)
		ds=0;
	lon0=(lon1+lon2)/2.0;
	zbar+=r*cos((lon1-lon0)*DEG_RAD)*cos(lat1*DEG_RAD)/4.0;
	xbar+=r*sin((lon1-lon0)*DEG_RAD)*cos(lat1*DEG_RAD)/4.0;
	zbar+=r*cos((lon1-lon0)*DEG_RAD)*cos(lat2*DEG_RAD)/4.0;
	xbar+=r*sin((lon1-lon0)*DEG_RAD)*cos(lat2*DEG_RAD)/4.0;
	zbar+=r*cos((lon2-lon0)*DEG_RAD)*cos(lat1*DEG_RAD)/4.0;
	xbar+=r*sin((lon2-lon0)*DEG_RAD)*cos(lat1*DEG_RAD)/4.0;
	zbar+=r*cos((lon2-lon0)*DEG_RAD)*cos(lat2*DEG_RAD)/4.0;
	xbar+=r*sin((lon2-lon0)*DEG_RAD)*cos(lat2*DEG_RAD)/4.0;
	ybar+=r*sin(lat1*DEG_RAD)/2.0;
	ybar+=r*sin(lat2*DEG_RAD)/2.0;
	lon0*=DEG_RAD;
	scale=6.2832*r*(lat2-lat1)/360.0/9.0;
	thick*=scale;
	drange=dmax-dmin;
	i0=(90.0-lat2)*12.0;
	i1=(90.0-lat1)*12.0;
	if(lon1<0.0)
		lon1+=360.0;
	j0=lon1*(double)12.0;
	if(lon2<0.0)
		lon2+=360.0;
	j1=lon2*(double)12.0;
	row=i1-i0+1;
	col=j1-j0+1;
	pts=(long)((row+samp-1)/samp)*(long)((col+samp-1)/samp);
	offset=(i0*dcol+(long)j0)*2;
	lat2=(double)90.0-(double)i0/12.0;
	lat1=(double)90.0-(double)i1/12.0;
	lon1=(double)j0/(double)12.0;
	lon2=(double)j1/(double)12.0;
	fpin=fopen_path("globe16.dat","rb");
	if(!fpin)
	{
		printf("Could not open 'globe16.dat' to read.\n\n");
		printf("If file is not on this disk include path, e.g.\n");
		printf("\netopo5 l:\\3d\\\n");
		exit(0);
	}
	sprintf(string,"%s.wvf",name);
	fpout=fopen(string,"wt");
	if(!fpout)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	printf("Do you want to reset scale and location parameters/ (y or n)\n\n");
	if(getch()=='y')
	{
		printf("X0 = %lf km\n",xbar);
		printf("Y0 = %lf km\n",ybar);
		printf("Z0 = %lf km\n",zbar);
		printf("scale = %lf km/in\n",scale);
		printf("lon0  = %lf deg\n",lon0/DEG_RAD);
		printf("Change X0?  (y or n)\n");
		if(getch()=='y')
		{
			printf("Give new X0\n");
			scanf("%lf",&xbar);
		}
		printf("Change Y0?  (y or n)\n");
		if(getch()=='y')
		{
			printf("Give new Y0\n");
			scanf("%lf",&ybar);
		}
		printf("Change Z0?  (y or n)\n");
		if(getch()=='y')
		{
			printf("Give new Z0\n");
			scanf("%lf",&zbar);
		}
		printf("Change Scale?  (y or n)\n");
		if(getch()=='y')
		{
			printf("Give new Scale\n");
			scanf("%lf",&scale);
		}
		printf("Change Front_Longitude?  (y or n)\n");
		if(getch()=='y')
		{
			printf("Give new Front_Longitude\n");
			scanf("%lf",&lon0);
			lon0*=DEG_RAD;
		}
	}
	fprintf(fpout,"# lat = %lf to %lf\n",lat1,lat2);
	fprintf(fpout,"# lon = %lf to %lf\n",lon1,lon2);
	fprintf(fpout,"# X0 = %lf km\n",xbar);
	fprintf(fpout,"# Y0 = %lf km\n",ybar);
	fprintf(fpout,"# Z0 = %lf km\n",zbar);
	fprintf(fpout,"# scale = %lf km/in\n",scale);
	fprintf(fpout,"# vertical exaggeration = %lf\n",zscale);
	fprintf(fpout,"# sample rate = %d\n",samp);
	fprintf(fpout,"# PROJECTION = spherical earth\n");
	fprintf(fpout,"# earth radius = %lf km\n",r);
	fprintf(fpout,"# %ld vertices\n",pts);
	if(ds!=0)
		fprintf(fpout,"# double sided figure -- thickness = %lf\n",thick);
	fprintf(fpout,"g %s\n",name);
	fprintf(fpout,"g\n");
	fseek(fpin,offset,SEEK_SET);
	printf("Writing point coordinatest to '%s.wvf'\n",name);
	for(i=0;i<row;i+=samp)
	{
		lat=(lat2-(double)i/12.0)*DEG_RAD;
		printf("%crow %d of %d  lat = %10.5lf",
				cr,i/samp+1,(row+samp-1)/samp,lat/DEG_RAD);
		fread((char *)Elev,sizeof(int),col,fpin);
		offset+=dcol*2*(long)(samp);
		fseek(fpin,offset,SEEK_SET);
		for(j=0;j<col;j+=samp)
		{
			lon=(lon1+(double)j/12.0)*DEG_RAD;
			rho=Elev[j];
			rho*=zscale;
			rho=r+rho/1000.0;
			x=(rho*sin(lon-lon0)*cos(lat)-xbar)/scale;
			y=(rho*sin(lat)-ybar)/scale;
			z=(rho*cos(lon-lon0)*cos(lat)-zbar)/scale;
			fprintf(fpout,"v %12.6lf %12.6lf %12.6lf\n",x,y,z);
			ppts+=1;
			if(ds!=0.0)
			{
				rho-=thick;
				x=(rho*sin(lon-lon0)*cos(lat)-xbar)/scale;
				y=(rho*sin(lat)-ybar)/scale;
				z=(rho*cos(lon-lon0)*cos(lat)-zbar)/scale;
				fprintf(fpout,"v %12.6lf %12.6lf %12.6lf\n",x,y,z);
				ppts+=1;
			}
			dval=Elev[j];
			dval-=dmin;
			dval*=255;
			E8[j]=dval/drange;
		}
	}
	lrow=row=(row+samp-1)/samp;
	lcol=col=(col+samp-1)/samp;
	mxpts=MAX_V;
	printf("\nWriting triangle point assignments to '%s.wvf'\n",name);
	if(ds==0)
	{		
		for(i=1;i<row;i++)
		{
			printf("%crow %d of %d",cr,i+1,row);
			ij[1][1]=(long)i*(long)col+2;
			ij[0][0]=(long)(i-1)*(long)col+1;
			ij[0][1]=ij[0][0]+1;
			ij[1][0]=ij[1][1]-1;
			for(j=1;j<col;j++)
			{
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][0],ij[0][1],ij[0][0]);
				fprintf(fpout,"f %ld %ld %ld\n",ij[1][0],ij[1][1],ij[0][1]);
				for(m=0;m<2;m++)
					for(n=0;n<2;n++)
						ij[m][n]+=1;
			}
		}
	}
	else
	{
		for(i=1;i<row;i++)
		{
			printf("%crow %d of %d",cr,i+1,row);
			lr=i;
			ij3[0][0][0]=lcol*l2*(lr-l1)+l1;	/* top upper left (1st cube) */
			ij3[0][1][0]=lcol*l2*(lr-l1)+l3;	/* top upper right (1st cube) */
			ij3[0][0][1]=lcol*l2*lr+l1;	/* top lower left (1st cube) */
			ij3[0][1][1]=lcol*l2*lr+l3;	/* top lower right (1st cube) */
	
			ij3[1][0][0]=lcol*l2*(lr-l1)+l2;	/* bottom upper left (1st cube) */
			ij3[1][1][0]=lcol*l2*(lr-l1)+l4;	/* bottom upper right (1st cube) */
			ij3[1][0][1]=lcol*l2*lr+l2;	/* bottom lower left (1st cube) */
			ij3[1][1][1]=lcol*l2*lr+l4;	/* bottom lower right (1st cube) */
			for(j=1;j<col;j++)
			{
				lc=j;
/* top surface */
			fprintf(fpout,"f %ld %ld %ld\n",ij3[0][0][1],ij3[0][1][0],ij3[0][0][0]);
			fprintf(fpout,"f %ld %ld %ld\n",ij3[0][0][1],ij3[0][1][1],ij3[0][1][0]);
/* bottom surface */
			fprintf(fpout,"f %ld %ld %ld\n",ij3[1][0][1],ij3[1][0][0],ij3[1][1][0]);
			fprintf(fpout,"f %ld %ld %ld\n",ij3[1][0][1],ij3[1][1][0],ij3[1][1][1]);
/* left edge */
				if(j==1)
				{
				fprintf(fpout,"f %ld %ld %ld\n",ij3[1][0][0],ij3[0][0][1],ij3[0][0][0]);
				fprintf(fpout,"f %ld %ld %ld\n",ij3[1][0][0],ij3[1][0][1],ij3[0][0][1]);
				}
/* right edge */
				if(j==col-1)
				{
				fprintf(fpout,"f %ld %ld %ld\n",ij3[1][1][1],ij3[0][1][0],ij3[0][1][1]);
				fprintf(fpout,"f %ld %ld %ld\n",ij3[1][1][1],ij3[1][1][0],ij3[0][1][0]);
				}
/* upper edge */
				if(i==1)
				{
				fprintf(fpout,"f %ld %ld %ld\n",ij3[1][1][0],ij3[1][0][0],ij3[0][0][0]);
				fprintf(fpout,"f %ld %ld %ld\n",ij3[1][1][0],ij3[0][0][0],ij3[0][1][0]);
				}
/* lower edge */
				if(i==row-1)
				{
				fprintf(fpout,"f %ld %ld %ld\n",ij3[1][0][1],ij3[0][1][1],ij3[0][0][1]);
				fprintf(fpout,"f %ld %ld %ld\n",ij3[1][0][1],ij3[1][1][1],ij3[0][1][1]);
				}
				for(l=0;l<2;l++)	/* move cube 1 to left */
					for(m=0;m<2;m++)
						for(n=0;n<2;n++)
							ij3[l][m][n]+=2;
			}
		}
	}
	fprintf(fpout,"# %ld elements\n",ppts);
	printf("\n");
	fclose(fpin);
	fclose(fpout);
	if(ppts>MAX_V)
	{
		printf(
"There are %ld verticies and 'WAVEPREP' can only handle %ld pts.\n\n",
			ppts,mxpts);

	}
	else
		execlp("waveprep.exe","waveprep",name,NULL);
}

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

int comb(char *in1,char *in2,char *out)

{
	int i,j,k;
	int row=1080,col=COL;
	char string[100];
	FILE *fplbl,*fpout,*fpin;

	sprintf(string,"%s16.lbl",out);
	fplbl=fopen(string,"wt");
	if(!fplbl)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	else
	{
		fprintf(fplbl,"FILE_TYPE          = IMAGE\n");
		fprintf(fplbl,"IMAGE_LINES        =  2160\n");
		fprintf(fplbl,"LINE_SAMPLES       =  4320\n");
		fprintf(fplbl,"SAMPLE_BITS        =  16\n");
		fprintf(fplbl,"IMAGE_POINTER      = '%s16.dat'\n",out);
		fclose(fplbl);
	}
	sprintf(string,"%s16.dat",out);
	fpout=fopen(string,"wb");
	if(!fpout)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	sprintf(string,"%s16.dat",in1);
	fpin=fopen_path(string,"rb");
	if(!fpin)
	{
		printf("Could not open '%s' to read.\n\n",string);
		printf("If file is not on this disk include path, e.g.\n");
		printf("\netopo5 l:\\3d\\\n");
		exit(0);
	}
	for(i=0;i<row;i++)
	{
		fread((char *)E16,sizeof(int),col,fpin);
		fwrite((char *)E16,sizeof(int),col,fpout);
	}
	fclose(fpin);
	sprintf(string,"%s16.dat",in2);
	fpin=fopen_path(string,"rb");
	if(!fpin)
	{
		printf("Could not open '%s' to read.\n\n",string);
		printf("If file is not on this disk include path, e.g.\n");
		printf("\netopo5 l:\\3d\\\n");
		exit(0);
	}
	for(i=0;i<row;i++)
	{
		fread((char *)E16,sizeof(int),col,fpin);
		fwrite((char *)E16,sizeof(int),col,fpout);
	}
	fclose(fpin);
	fclose(fpout);

	sprintf(string,"%s8.lbl",out);
	fplbl=fopen(string,"wt");
	if(!fplbl)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	else
	{
		fprintf(fplbl,"FILE_TYPE          = IMAGE\n");
		fprintf(fplbl,"IMAGE_LINES        =  2160\n");
		fprintf(fplbl,"LINE_SAMPLES       =  4320\n");
		fprintf(fplbl,"IMAGE_POINTER      = '%s8.dat'\n",out);
		fclose(fplbl);
	}
	sprintf(string,"%s8.dat",out);
	fpout=fopen(string,"wb");
	if(!fpout)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	sprintf(string,"%s8.dat",in1);
	fpin=fopen_path(string,"rb");
	if(!fpin)
	{
		printf("Could not open '%s' to read.\n\n",string);
		printf("If file is not on this disk include path, e.g.\n");
		printf("\netopo5 l:\\3d\\\n");
		exit(0);
	}
	for(i=0;i<row;i++)
	{
		fread((char *)E8,sizeof(char),col,fpin);
		fwrite((char *)E8,sizeof(char),col,fpout);
	}
	fclose(fpin);
	sprintf(string,"%s8.dat",in2);
	fpin=fopen_path(string,"rb");
	if(!fpin)
	{
		printf("Could not open '%s' to read.\n\n",string);
		printf("If file is not on this disk include path, e.g.\n");
		printf("\netopo5 l:\\3d\\\n");
		exit(0);
	}
	for(i=0;i<row;i++)
	{
		fread((char *)E8,sizeof(char),col,fpin);
		fwrite((char *)E8,sizeof(char),col,fpout);
	}
	fclose(fpin);
	fclose(fpout);

}



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

int read_ascii(char *name)

{
	int i,j,k,l,m,n;
	FILE *fpin,*fpout16,*fpout8,*fplbl;
	int end,max=-32000,min=32000;
	char cr=0xd,string[100];
	long dmax=9000,dmin=-12000,drange,dval,dnum=0;
	char *cb;

	drange=dmax-dmin;
	sprintf(string,"%s.dat",name);
	fpin=fopen_path(string,"rt");
	if(!fpin)
	{
		printf("Could not open '%s' to read.\n\n",string);
		printf("If file is not on this disk include path, e.g.\n");
		printf("\netopo5 l:\\3d\\\n");
		exit(0);
	}
	sprintf(string,"%s16.lbl",name);
	fplbl=fopen(string,"wt");
	if(!fplbl)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	else
	{
		fprintf(fplbl,"FILE_TYPE          = IMAGE\n");
		fprintf(fplbl,"IMAGE_LINES        =  1080\n");
		fprintf(fplbl,"LINE_SAMPLES       =  4320\n");
		fprintf(fplbl,"SAMPLE_BITS        =  16\n");
		fprintf(fplbl,"IMAGE_POINTER      = '%s16.dat'\n",name);
		fclose(fplbl);
	}
	sprintf(string,"%s16.dat",name);
	fpout16=fopen(string,"wb");
	if(!fpout16)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	sprintf(string,"%s8.lbl",name);
	fplbl=fopen(string,"wt");
	if(!fplbl)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	else
	{
		fprintf(fplbl,"FILE_TYPE          = IMAGE\n");
		fprintf(fplbl,"IMAGE_LINES        =  1080\n");
		fprintf(fplbl,"LINE_SAMPLES       =  4320\n");
		fprintf(fplbl,"IMAGE_POINTER      = '%s8.dat'\n",name);
		fclose(fplbl);
	}
	sprintf(string,"%s8.dat",name);
	fpout8=fopen(string,"wb");
	if(!fpout8)
	{
		printf("Could not open '%s' to write.\n\n",string);
		printf(
"Transfer to write disk and use path on command line,e.g.\n\n");
		printf("etopo5 l:\\3d\\\n\n");
		exit(0);
	}
	end=0;
	l=0;
	do
	{
		for(k=0;k<COL;k++)
		{
			if(fscanf(fpin,"%d",Elev+k)<1)
				end=1;
			if(end==0)
			{
				sscanf(cb,"%d",Elev+k);
				dnum+=1;
				if(Elev[k]<min)
					min=Elev[k];
				if(Elev[k]>max)
					max=Elev[k];
				dval=Elev[k];
				dval-=dmin;
				dval*=255;
				E8[k]=dval/drange;
			}
		}
		fwrite((char *)Elev,sizeof(int),COL,fpout16);
		fwrite((char *)E8,sizeof(char),COL,fpout8);
		printf("%c%6d %6d %6d",cr,++l,min,max);
	}while(end==0);
	printf("\n\n%s max = %d  min = %d\n\n",name,max,min);
	printf("number of points = %ld\n",dnum);
}