#include <stdio.h>
#include <CoreFoundation/CoreFoundation.h>

int readGLAS(Float32 *buf, Float32 scaleFactor, const char *fname);

int main (int argc, const char * argv[]) {
	
	const int nx = 11352;	// These are GLAS antarctica grid dimensions
	const int ny =  9368;
	const int nn = nx * ny;
	
	Float32 *el, *lat, *lon;
	
	el  = malloc(sizeof(Float32)*nn);
	lat = malloc(sizeof(Float32)*nn);
	lon = malloc(sizeof(Float32)*nn);
	
	// Get 4 byte fixed point representation from file
	// Convert from 4 byte fixed point to 4 byte float using scale factor
	
	readGLAS( el, 1.e+3, "/Users/richard/SRTM15/icesat/NSIDC_Ant500m_egm96_elev_cm.dat");
	readGLAS(lat, 1.e+6, "/Users/jjbecker/Desktop/ICESAT_DEM/glas_ant_500m_Lat_udegN_9368x11352.dat");
	readGLAS(lon, 1.e+6, "/Users/jjbecker/Desktop/ICESAT_DEM/glas_ant_500m_Lon_udegE_9368x11352.dat");
	
	// Write out 4 byte floating point representation
	
	int fd = creat("glasFloat32XYZ.dat", 0666);
	int cnt = nn;
	//cnt = ny; // debug
	while (cnt--) {
		// Don't write NaN as "xyz2grd" will do it for us...
		if (*el > -999.) {
			if (*lon >= 180.) {
				*lon -= 360.;
			}
 			printf(" %f %f %f\n", *lat, *lon, *el);	// debug
			write(fd, lon++, 4);
			write(fd, lat++, 4);
			write(fd,  el++, 4);
		} else {
 			printf(" %f %f %f SKIPPED\n", *lat, *lon, *el);	// debug
			lat++; lon++; el++;
		}
	}
	close(fd);
	return 0;
}

int readGLAS(Float32 *buf, Float32 scaleFactor, const char *fname) {
	int fd = open(fname, 0);
	int cnt = 0;
	int numBytesRead = 0;
	int tmp;
	
	if (fd < 0) {
		printf("\n%s not found; exiting\n\n",fname);
		exit(1);
	}
	/*
	 This GLAS data file was written by FORTRAN. 
	 That's right: FORTRAN.
	 Unbelievable?
	 So the data is column major (i.e. all of col(1), then col(2) and so on)
	 */
	
	while ((numBytesRead = read(fd, &tmp, 4)) == 4) {
		tmp = (int)CFSwapInt32BigToHost(tmp);
		*buf++ = (float)tmp/scaleFactor;
		cnt++;
		//if (cnt > 9368) return cnt; // debug
	}
	return cnt;
}