/*
	elev.c


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

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

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

main()
{
	int i,j,k;
	int row,col;
	char ans,string[100],cr=0xd;
	double lat1,lat2,lon1,lon2,clat,r,km_in,dlon,dlat;
	double lt1,ln1,lt2,ln2,thick=0.5;
	long max,min;

	do
	{
		printf("   1 -- read North ASCII file and make 16-bit '.img'\n");
		printf("   2 -- read South ASCII file and make 16-bit '.img'\n");
		printf("   3 -- combine N and S images\n");
		printf("   4 -- make a sub image in Wave front format (round)\n");
		printf("   5 -- make a sub image in Wave front format (flat/solid)\n");
		printf("   6 -- make a sub image in Wave front format (flat/skin)\n");
		printf("   7 -- make a sub image set in Wave front format (flat/solid)\n");

		printf("\n   a -- check max min of set\n");

		printf("\n   x -- exit program\n\n");

		ans=getch();
		if(ans=='1')
			read_ascii("elev_n");
		if(ans=='2')
			read_ascii("elev_s");
		if(ans=='3')
			comb("elev_n","elev_s","globe");
		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 earth radius in kilometers (640 gives 10 x relief)\n");
			scanf("%lf",&r);
			wf(lat1,lat2,lon1,lon2,r);
		}
		if(ans=='5')
		{
			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=='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 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=='7')
		{
			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=='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');
}

/**********************************************************************
**
**	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];
	double scale=cos(clat*DEG_RAD);
	long nm,cl,rw,lr,lc,lrow,lcol,l2=2,l1=1,l3=3,l4=4;

	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);
		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);
		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;
	fpin=fopen("globe16.dat","rb");
	if(!fpin)
	{
		printf("Could not open 'globe16.dat' to read.\n\n");
		exit(0);
	}
	fpout=fopen(name,"wt");
	if(!fpout)
	{
		printf("Could not open '%s' to write.\n\n",name);
		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++)
	{
		lat=(lat2-(double)i/12.0);
		y=-(lat2-lat)/(double)360.0*E_CIR*km_in-4.5;
		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-lon1)/(double)360.0*E_CIR*scale*km_in-4.5;
			z=rho*r*km_in;
			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);
	}
	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;

	*max=-1000000000,*min=1000000000;
	i0=(90.0-lat2)*12.0;
	if(i0<0||i0>=COL)
	{
		printf("maximum latitude is not valid.\n\n");
		return(-1);
	}
	i1=(90.0-lat1)*12.0;
	if(i1<0||i1>=COL)
	{
		printf("minimum latitude is not valid.\n\n");
		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");
		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");
		return(-1);
	}
	row=i1-i0+1;
	col=j1-j0+1;
	if(row<0||row>ROW)
	{
		printf("Latitude range is not valid.\n\n");
		printf(
			"You may have crossed the 0 deg meridian which is not allowed.\n\n");
		return(-1);
	}
	if(col<0||col>COL)
	{
		printf("longitude range is not valid.\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("globe16.dat","rb");
	if(!fpin)
	{
		printf("Could not open 'globe16.dat' to read.\n\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-4.5;
		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);
}




/**********************************************************************
**
**	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;
	long offset,pts,ij[2][2];
	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;
	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);
		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);
		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("globe16.dat","rb");
	if(!fpin)
	{
		printf("Could not open 'globe16.dat' to read.\n\n");
		exit(0);
	}
	fpout=fopen(name,"wt");
	if(!fpout)
	{
		printf("Could not open '%s' to write.\n\n",name);
		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,"# %7.3lf km/in\n",(double)1.0/km_in);
	fprintf(fpout,"# Z-Range  %ld to %ld\n",dmin,dmax);
	fprintf(fpout,"# Z-Mid = 0\n");
	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=-(lat2-lat)/(double)360.0*E_CIR*km_in-4.5;
		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-lon1)/(double)360.0*E_CIR*scale*km_in-4.5;
			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-4.5;
	for(j=0;j<cl;j++)
	{
		lon=(lon1+(double)j/12.0);
		x=(lon-lon1)/(double)360.0*E_CIR*scale*km_in-4.5;
		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-4.5;
	for(i=0;i<rw;i++)
	{
		lat=(lat2-(double)i/12.0);
		y=-(lat2-lat)/(double)360.0*E_CIR*km_in-4.5;
		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-4.5;
	for(j=cl;j>0;j--)
	{
		lon=(lon1+(double)j/12.0);
		x=(lon-lon1)/(double)360.0*E_CIR*scale*km_in-4.5;
		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-4.5;
	for(i=rw;i>0;i--)
	{
		lat=(lat2-(double)i/12.0);
		y=-(lat2-lat)/(double)360.0*E_CIR*km_in-4.5;
		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)

{
	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,*fpout,*fplbl,*fpdat;
	long dmax=9000,dmin=-12000,drange,dval,dnum=0;
	char string[100];

	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);
		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);
		exit(0);
	}

	pts=(long)row*(long)col;
	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("globe16.dat","rb");
	if(!fpin)
	{
		printf("Could not open 'globe16.dat' to read.\n\n");
		exit(0);
	}
	fpout=fopen("globe.wvf","wt");
	if(!fpout)
	{
		printf("Could not open 'globe.wvf' to write.\n\n");
		exit(0);
	}
	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");
	fprintf(fpout,"g\n");
	fseek(fpin,offset,SEEK_SET);
	for(i=0;i<row;i++)
	{
		lat=(lat2-(double)i/12.0)*DEG_RAD;
		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)*DEG_RAD;
			rho=Elev[j];
			rho=r+rho/1000.0;
			x=rho*cos(lon)*cos(lat)-4.5;
			y=rho*sin(lon)*cos(lat)-4.5;
			z=rho*sin(lat);
			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);
	}
	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;
		}
	}
	pts=(long)(row-1)*(long)(col-1)*2;
	fprintf(fpout,"# %ld elements\n",pts);
	fclose(fpin);
	fclose(fpout);
	fclose(fpdat);
}

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

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);
		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);
		exit(0);
	}
	sprintf(string,"%s16.dat",in1);
	fpin=fopen(string,"rb");
	if(!fpin)
	{
		printf("Could not open '%s' to read.\n\n",string);
		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(string,"rb");
	if(!fpin)
	{
		printf("Could not open '%s' to read.\n\n",string);
		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);
		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);
		exit(0);
	}
	sprintf(string,"%s8.dat",in1);
	fpin=fopen(string,"rb");
	if(!fpin)
	{
		printf("Could not open '%s' to read.\n\n",string);
		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(string,"rb");
	if(!fpin)
	{
		printf("Could not open '%s' to read.\n\n",string);
		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(string,"rt");
	if(!fpin)
	{
		printf("Could not open '%s' to read.\n\n",string);
		exit(0);
	}
	sprintf(string,"%s16.lbl",name);
	fplbl=fopen(string,"wt");
	if(!fplbl)
	{
		printf("Could not open '%s' to write.\n\n",string);
		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);
		exit(0);
	}
	sprintf(string,"%s8.lbl",name);
	fplbl=fopen(string,"wt");
	if(!fplbl)
	{
		printf("Could not open '%s' to write.\n\n",string);
		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);
		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);
}