/* -*- MaTX -*-
 *
 *  NAME
 *      {M1,M2,,...,Mn,Names} = matlab_read(filename)
 *        Matrix M1,M2,...,Mn;
 *        List Names = {name1, name2, ..., namen};
 *        String name1, name2, ..., namen;
 *        String filename;
 *
 *      matlab_write(filename, matrices)
 *        String filename;
 *        List matrices;
 *
 *      matlab_write(filename, matrices, type)
 *        String filename;
 *        List matrices;
 *        String type;
 *
 *  DESCRIPTION
 *      matlab_read() load matrices from a file saved as matlab 
 *      mat-format (version 4)
 *      matlab_write() save matrices to a file as matlab mat-format
 *      (version 4)
 *
 *  SEE ALSO
 *      fread and fwrite
 */

Func void matlab_write(filename, matrices, type, ...)
	String filename;
	List matrices;
	String type;
{
	Integer i, fd;

	void matlab_write_data();

	error(nargchk(2, 3, nargs, "matlab_write"));

	if (nargs == 2) {
		type = "double";
	}

	if ((fd = fopen(filename, "bw")) < 0) {
		error("matlab_write(): Can't open %s\n", filename);
	}

	for (i = 1; i <= length(matrices); i++) {
		matlab_write_data(fd, matrices(i,Matrix), sprintf("mat%d",i), type);
	}

	fclose(fd);
}

Func void matlab_write_data(fd, a, name, prec)
	Integer fd;
	Matrix a;
	String name;
	String prec;
{
	Integer type, rows, cols, imagf, namelen;
	Integer M,O,P,T;
	Matrix a2;
	
	M = machine_endian();
	O = 0;
	if (prec == "double") {
		P = 0;
	} else if (prec == "float") {
		P = 1;
	} else if (prec == "long") {
		P = 2;
	} else if (prec == "short") {
		P = 3;
	} else if (prec == "unsigned short") {
		P = 4;
	} else if (prec == "unsigned char") {
		P = 5;
	} else {
		error("matlab_write_data(): Incorrect data precision\n");
	}
	T = 0;

	type = M*1000 + O*100 + P*10 + T;
	rows = Rows(a);
	cols = Cols(a);
	imagf = iscomplex(a);
	namelen = length(name) + 1;

	// write header
	fwrite(fd, [type rows cols imagf namelen], "unsigned long");
	fwrite(fd, [abs(name), [0]], "unsigned char");

	a2 = reshape(trans(a), rows, cols);
	if (imagf) {
		fwrite(fd, Re(a2), prec);
		fwrite(fd, Im(a2), prec);
	} else {
		fwrite(fd, a2, prec);
	}
}

Func List matlab_read(filename)
	String filename;
{
	Integer fd, flip;
	Matrix a;
	List matrices, names;
	String mname;

	Integer check_data_endian();
	List matlab_read_data();

	flip = check_data_endian(filename);

	if ((fd = fopen(filename, "br")) < 0) {
		error("matlab_read(): Can't open %s\n", filename);
	}
	
	matrices = {};
	names = {};
	while (1) {
		{a,mname} = matlab_read_data(fd, flip);
		if (length(a) == 0) {
			break;
		}
		appendlist(matrices,a);
		appendlist(names,mname);
	}

	fclose(fd);
	appendlist(matrices,names);
	return matrices;
}

Func List matlab_read_data(fd, flip)
	Integer fd, flip;
{
	Integer type, rows, cols, imagf, namelen;
	Integer M,O,P,T;
	Matrix header, data, rdata, idata;
	String mname;

	// Read header
	header = fread(fd, 5, "unsigned long", flip);
	if (length(header) == 0) {
		return {[],""};
	}

	// Header
	type = Integer(header(1));
	rows = Integer(header(2));
	cols = Integer(header(3));
	imagf = Integer(header(4));
	namelen = Integer(header(5));

	M = type/1000;
	O = rem(type, 1000)/100;
	P = rem(type, 100)/10;
	T = rem(type, 10);
	
	if (T == 1) {
		error("matlab_read_data(): Can't read text matrix\n");
	} else if (T == 2) {
		error("matlab_read_data(): Can't read sparce matrix\n");
	}

	mname = String(fread(fd, namelen, "char"));
	if (length(mname) == 0) {
		error("matlab_read_data(): Can't read name\n");
	}
	
	if (imagf) {
		imagf = 2;
	} else {
		imagf = 1;
	}

	switch (P) {
	  case 0:
		data = fread(fd, imagf*rows*cols, "double", flip);
		break;
	  case 1:
		data = fread(fd, imagf*rows*cols, "float", flip);
		break;
	  case 2:
		data = fread(fd, imagf*rows*cols, "long", flip);
		break;
	  case 3:
		data = fread(fd, imagf*rows*cols, "short", flip);
		break;
	  case 4:
		data = fread(fd, imagf*rows*cols, "unsigned short", flip);
		break;
	  case 5:
		data = fread(fd, imagf*rows*cols, "unsigned char", flip);
		break;
	  default:
		error("matlab_read_data(): Incorrect precision\n");
		break;
	}
	
	if (length(data) == 0) {
		error("matlab_read_data(): Too few data\n");
	}
	
	if (imagf == 1) {
		data = trans(reshape(data, cols, rows));
	} else {
		rdata = data(1:rows*cols);
		idata = data(rows*cols+1:2*rows*cols);
		data = (rdata,idata);
		data = trans(reshape(data, cols, rows));
	}

	return {data, mname};
}

Func Integer check_data_endian(filename)
	String filename;
{
	Integer fd, flip;
	Matrix header;

	flip = 0;

	if ((fd = fopen(filename, "br")) < 0) {
		error("check_data_endian(): Can't open %s\n", filename);
	}

	// Read header
	header = fread(fd, 5, "unsigned long", flip);
	if (length(header) == 0) {
		error("check_data_endian(): Can't read header\n");
	}
	
	if (header(1) < 0 || 9999 < header(1) ||
		(header(1) == 0 && machine_endian() != 0)) {
		// Machine endian is dirrent
		flip = 1;
	}

	fclose(fd);
	return flip;
}