/*
 Simple Data Format
 http://solarmuri.ssl.berkeley.edu/~fisher/public/software/SDF
 Copyright (C) 2006 University of California

 This is free software; you can redistribute it and/or
 modify it under the terms of the GNU Lesser General Public
 License as published by the Free Software Foundation;
 either version 2.1 of the License, or (at your option) any later version.

 This software is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 See the GNU Lesser General Public License for more details.

 To view the GNU Lesser General Public License visit
 http://www.gnu.org/copyleft/lesser.html
 or write to the Free Software Foundation, Inc.,
 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <sdf_subs.h>
void sdf_write_f77__(char *fnamef77, char *labf77, char *dtf77, intf *nbpwf77, 
       intf *ndimf77, intf *dimsf77, void *dataf77,  intf len_fnamef77, 
       intf len_labf77, intf len_dtf77, intf len_dataf77 )
{
sdf_write_f77( fnamef77, labf77, dtf77, nbpwf77, 
      ndimf77, dimsf77, dataf77,  len_fnamef77, 
       len_labf77, len_dtf77, len_dataf77 );
return;
}
void sdf_write_f77_(char *fnamef77, char *labf77, char *dtf77, intf *nbpwf77, 
       intf *ndimf77, intf *dimsf77, void *dataf77,  intf len_fnamef77, 
       intf len_labf77, intf len_dtf77, intf len_dataf77 )
{
sdf_write_f77( fnamef77, labf77, dtf77, nbpwf77, 
      ndimf77, dimsf77, dataf77,  len_fnamef77, 
       len_labf77, len_dtf77, len_dataf77 );
return;
}
void sdf_write_f77(char *fnamef77, char *labf77, char *dtf77, intf *nbpwf77, 
       intf *ndimf77, intf *dimsf77, void *dataf77,  intf len_fnamef77, 
       intf len_labf77, intf len_dtf77, intf len_dataf77 )
{

/*
NOTE -- Fortran character
strings are not passed as a simple pointer, as one might naively guess.  
There's a 2nd hidden integer argument related to the character string length.
Think maybe I got this problem solved now....
*/

data_id *id;
i4 nlab,i;
char *fname;
char *datatype;
char *label;

/* 1st, print out the sizes of f77 char arrays (debug) */

/*
printf("len_fnamef77 = %d\n",len_fnamef77);
printf("len_labf77 = %d\n",len_labf77);
printf("len_dtf77 = %d\n",len_dtf77);
printf("len_dataf77 = %d\n",len_dataf77);
*/

/* Convert fortran strings fnamef77,labf77,dtf77 to C-strings fname,label,
   datatype */

fs2cs(&fname,fnamef77,len_fnamef77);
fs2cs(&label,labf77,len_labf77);
fs2cs(&datatype,dtf77,len_dtf77);

/* Debug - print out C-strings fname,label and character variable datatype */

/*
printf("fname = %s\n",fname);
printf("label = %s\n",label);
printf("datatype = %c\n",*datatype);
*/

id=malloc(sizeof(data_id));
nlab=strlen(label);
id->label=calloc(nlab+1,1);

for (i=0;i<nlab;i++)
{
	*(id->label+i)=*(label+i);
}

id->nbpw=(i4) *nbpwf77;
id->datatype=*datatype;
id->order=0; /* value shouldn't matter */
id->ndim=(i4) *ndimf77;
id->dims=malloc(id->ndim * sizeof(pos));

/* I wonder if the f77 integer values of the dimensions shouldn't be
   assumed to be of regular integer type rather than (pos).  
   They could be recast as (pos) inside this fn.  OK, now implemented. */

for (i=0;i<id->ndim;i++)
{
	*(id->dims+i)= (pos) *(dimsf77+id->ndim-1-i); /* flip dims for f77 */
}
sdf_write(fname, id, dataf77);
free(label);
free(fname);
free(datatype);
free(id->label);
free(id->dims);
free(id);
}

void sdf_read_f77__(char *fnamef77, intf *iorder_f77, char *labf77, 
     char *dtf77, intf *nbpwf77, intf *ndimf77, intf *dimsf77, void *dataf77, 
     intf len_fnamef77, intf len_labf77, intf len_dtf77, intf len_dataf77)
{
sdf_read_f77(fnamef77, iorder_f77, labf77, dtf77, 
     nbpwf77, ndimf77, dimsf77, dataf77, 
     len_fnamef77, len_labf77, len_dtf77, len_dataf77);
return;
}

void sdf_read_f77_(char *fnamef77, intf *iorder_f77, char *labf77, char *dtf77, 
     intf *nbpwf77, intf *ndimf77, intf *dimsf77, void *dataf77, 
     intf len_fnamef77, intf len_labf77, intf len_dtf77, intf len_dataf77)
{
sdf_read_f77(fnamef77, iorder_f77, labf77, dtf77, 
     nbpwf77, ndimf77, dimsf77, dataf77, 
     len_fnamef77, len_labf77, len_dtf77, len_dataf77);
return;
}

void sdf_read_f77(char *fnamef77, intf *iorder_f77, char *labf77, char *dtf77, 
     intf *nbpwf77, intf *ndimf77, intf *dimsf77, void *dataf77, 
     intf len_fnamef77, intf len_labf77, intf len_dtf77, intf len_dataf77)

{
data_id *id;
i4 nlab;
char *fname;
FILE *fp;
pos datasize, curpos;
char *sdfstr="SDF format";
char temp[301], testid[300];
i4 i,ise,ibe, iorder, order, ndatasets;
pos hdrsize;

/* Make certain (pos) is defined as 64bit int, otherwise, exit */
test_sizes(); 

/* DEBUG - print out lengths of fortran stings (hidden args) */
/*
printf("len_fnamef77 = %d\n",len_fnamef77);
printf("len_labf77 = %d\n",len_labf77);
printf("len_dtf77 = %d\n",len_dtf77);
printf("len_dataf77 = %d\n",len_dataf77);
*/

/* Convert fnamef77 to C-string fname */
fs2cs(&fname,fnamef77,len_fnamef77);
order = (i4) *iorder_f77;
/* DEBUG - print out C-string version of filename fname */
/*
printf("fname = %s\n",fname);
*/

ibe = is_big_endian ();
ise = 0;
if (ibe == 0) ise = 1;    /* set flag for byteswapping if small endian */

/* Initially, open the file as read-only */

fp=fopen(fname,"rb");

if (fp == NULL) /* does file fname exist? */
{
	printf("sdf_read_f77: file %s does not exist, exiting\n",fname);
	fflush(stdout);
	exit(1);
}
else
{
	fgets(testid,strlen(sdfstr)+1,fp);
	/* Now, test to make sure the SDF file identifier string is there */
	if(strncmp(sdfstr,testid,strlen(sdfstr)))
	{
		fclose(fp);
		printf("sdf_read_f77: File %s is not an SDF file, exiting\n",
			fname);
		fflush(stdout);
		exit(1);
	}
}
/* read the number of datasets in the file: */

fseeko(fp, strlen(sdfstr)+1 + 2*sizeof(pos) , SEEK_SET);
fread(&ndatasets,sizeof(i4),1,fp);
if(ise) byteswap((void *)&ndatasets,(pos)1,sizeof(i4));
/*  This print statement possibly important, but nuke for now */

/*
printf("sdf_read_f77: existing datasets = %d\n",ndatasets);
fflush(stdout);
*/

/*
printf("sdf_read_f77: desired dataset order = %d\n",order);
*/
if(order > (ndatasets-1))
{
	printf("sdf_read_f77:  Desired dataset order %d > ndatasets-1 = %d\n",
             order, ndatasets-1);
	fflush(stdout);
	exit(1);
}
/* read in header size */
fread(&hdrsize,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrsize,(pos)1,sizeof(pos));
/* debug */
/*
output_int64(temp,hdrsize);
printf("sdf_read_f77: hdrsize = %s\n",temp);
fflush(stdout);
*/

/* read the string data of all the datasets, print it out */
curpos = (pos) strlen(sdfstr)+1+sizeof(pos)+hdrsize;
/*
        printf("DEBUG - sdf_read_f77: curpos = %qu\n",curpos);
*/
for (iorder=0; iorder <= order; iorder++)
{
/*
        fgets(temp,300,fp);
        printf("iorder = %d , header string = %s",iorder,temp);
*/

        id=string_to_data(fp);
	datasize=data_size(id);


        if(iorder < order)
	{
		if(id->datatype != 'c') /* test for complex var */
		{
			curpos += (pos) datasize*(id->nbpw);
		}
		else
		{
			curpos += (pos) datasize*2*(id->nbpw);
			/* kludge for complex variables */
		}
		free(id->label);
		free(id->dims);
		free(id);
	}
	else
	{

	output_int64(temp,datasize);

        /*
	 These print statements possibly important, but nuke for now
        printf("sdf_read_f77: dataset %d, size %selements, descriptor: ",
               iorder,temp);
	fflush(stdout);
        data_to_string(id,stdout);
	printf("\n");
	fflush(stdout);
        */
	}
}

fseeko(fp,curpos,SEEK_SET);                /* reposition file pointer  */

if (id->datatype != 'c') /* test for complex var */
{
	fread(dataf77, id->nbpw , datasize,fp);      /* read the data */
	if(ise) byteswap(dataf77,datasize,id->nbpw); /* byteswap if neeeded */
}
else
{
	/* kludge for reading complex variable data */
	fread(dataf77, id->nbpw , 2*datasize,fp);      /* read the data */
	if(ise) byteswap(dataf77,2*datasize,id->nbpw); /* byteswap if neeeded */
}

nlab=strlen(id->label);

/* now convert id->label into fortran string labf77 */
cs2fs(id->label,labf77,len_labf77);

*dtf77 = id->datatype;
*nbpwf77 = (intf) id->nbpw;
*ndimf77 = (intf) id->ndim;
for (i=0; i < id->ndim; i++)
{
	/* flip order for f77 */
	*(dimsf77+i) = (intf) *(id->dims+id->ndim-1 -i); 
}
datasize=data_size(id);

free(fname);
free(id->label);
free(id->dims);
free(id);
fclose(fp);

}

void sdf_delete_f77__(char *fnamef77, intf* idelete, intf len_fnamef77)
{
sdf_delete_f77(fnamef77, idelete, len_fnamef77);
return;
}
void sdf_delete_f77_(char *fnamef77, intf* idelete, intf len_fnamef77)
{
sdf_delete_f77(fnamef77, idelete, len_fnamef77);
return;
}
void sdf_delete_f77(char *fnamef77, intf* idelete, intf len_fnamef77)
{

/*
NOTE -- Fortran character
strings are not passed as a simple pointer, as one might naively guess.  
There's a 2nd hidden integer argument related to the character string length.
Think maybe I got this problem solved now....
*/

char *fname;

/* 1st, print out the sizes of f77 char arrays (debug) */

/*
printf("len_fnamef77 = %d\n",len_fnamef77);
printf("len_labf77 = %d\n",len_labf77);
printf("len_dtf77 = %d\n",len_dtf77);
printf("len_dataf77 = %d\n",len_dataf77);
*/

/* Convert fortran strings fnamef77,labf77,dtf77 to C-strings fname,label,
   datatype */

fs2cs(&fname,fnamef77,len_fnamef77);

/* Debug - print out C-strings fname,label and character variable datatype */

/*
printf("fname = %s\n",fname);
printf("label = %s\n",label);
printf("datatype = %c\n",*datatype);
*/


/* I wonder if the f77 integer values of the dimensions shouldn't be
   assumed to be of regular integer type rather than (pos).  
   They could be recast as (pos) inside this fn.  OK, now implemented. */

sdf_delete(fname, (i4) *idelete);
free(fname);
}

void sdf_insert_f77__(char *fnamef77, intf* insert, char *labf77, 
       char *dtf77, intf *nbpwf77, intf *ndimf77, intf *dimsf77, 
       void *dataf77,  intf len_fnamef77, intf len_labf77, intf len_dtf77, 
       intf len_dataf77 )
{
sdf_insert_f77(fnamef77, insert, labf77, 
       dtf77, nbpwf77, ndimf77, dimsf77, 
       dataf77,  len_fnamef77, len_labf77, len_dtf77, 
       len_dataf77 );
return;
}
void sdf_insert_f77_(char *fnamef77, intf* insert, char *labf77, 
       char *dtf77, intf *nbpwf77, intf *ndimf77, intf *dimsf77, 
       void *dataf77,  intf len_fnamef77, intf len_labf77, intf len_dtf77, 
       intf len_dataf77 )
{
sdf_insert_f77(fnamef77, insert, labf77, 
       dtf77, nbpwf77, ndimf77, dimsf77, 
       dataf77,  len_fnamef77, len_labf77, len_dtf77, 
       len_dataf77 );
return;
}

void sdf_insert_f77(char *fnamef77, intf* insert, char *labf77, 
       char *dtf77, intf *nbpwf77, intf *ndimf77, intf *dimsf77, 
       void *dataf77,  intf len_fnamef77, intf len_labf77, intf len_dtf77, 
       intf len_dataf77 )
{

/*
NOTE -- Fortran character
strings are not passed as a simple pointer, as one might naively guess.  
There's a 2nd hidden integer argument related to the character string length.
Think maybe I got this problem solved now....
*/

data_id *id_new;
i4 nlab,i;
char *fname;
char *datatype;
char *label;

/* 1st, print out the sizes of f77 char arrays (debug) */

/*
printf("len_fnamef77 = %d\n",len_fnamef77);
printf("len_labf77 = %d\n",len_labf77);
printf("len_dtf77 = %d\n",len_dtf77);
printf("len_dataf77 = %d\n",len_dataf77);
*/

/* Convert fortran strings fnamef77,labf77,dtf77 to C-strings fname,label,
   datatype */

fs2cs(&fname,fnamef77,len_fnamef77);
fs2cs(&label,labf77,len_labf77);
fs2cs(&datatype,dtf77,len_dtf77);

/* Debug - print out C-strings fname,label and character variable datatype */

/*
printf("fname = %s\n",fname);
printf("label = %s\n",label);
printf("datatype = %c\n",*datatype);
*/

id_new=malloc(sizeof(data_id));
nlab=strlen(label);
id_new->label=calloc(nlab+1,1);

for (i=0;i<nlab;i++)
{
	*(id_new->label+i)=*(label+i);
}

id_new->nbpw=(i4) *nbpwf77;
id_new->datatype=*datatype;
id_new->order=0; /* value shouldn't matter */
id_new->ndim=(i4) *ndimf77;
id_new->dims=malloc(id_new->ndim * sizeof(pos));

/* I wonder if the f77 integer values of the dimensions shouldn't be
   assumed to be of regular integer type rather than (pos).  
   They could be recast as (pos) inside this fn.  OK, now implemented. */

for (i=0;i<id_new->ndim;i++)
{
	*(id_new->dims+i)= (pos) *(dimsf77+id_new->ndim-1-i); /*flip dims f77*/
}
sdf_insert(fname, (i4) *insert, id_new, dataf77);
free(label);
free(fname);
free(datatype);
free(id_new->label);
free(id_new->dims);
free(id_new);
}

void sdf_replace_f77__(char *fnamef77, intf * replace, char *labf77, 
       char *dtf77, intf  *nbpwf77, intf  *ndimf77, intf  *dimsf77, 
       void *dataf77,  intf  len_fnamef77, intf  len_labf77, intf  
       len_dtf77, intf  len_dataf77 )
{
sdf_replace_f77(fnamef77, replace, labf77, 
       dtf77, nbpwf77, ndimf77, dimsf77, 
       dataf77, len_fnamef77, len_labf77, len_dtf77, len_dataf77 );
return;
}
void sdf_replace_f77_(char *fnamef77, intf * replace, char *labf77, 
       char *dtf77, intf  *nbpwf77, intf  *ndimf77, intf  *dimsf77, 
       void *dataf77,  intf  len_fnamef77, intf  len_labf77, intf  
       len_dtf77, intf  len_dataf77 )
{
sdf_replace_f77(fnamef77, replace, labf77, 
       dtf77, nbpwf77, ndimf77, dimsf77, 
       dataf77, len_fnamef77, len_labf77, len_dtf77, len_dataf77 );
return;
}
void sdf_replace_f77(char *fnamef77, intf * replace, char *labf77, 
       char *dtf77, intf  *nbpwf77, intf  *ndimf77, intf  *dimsf77, 
       void *dataf77,  intf  len_fnamef77, intf  len_labf77, intf  
       len_dtf77, intf  len_dataf77 )
{

/*
NOTE -- Fortran character
strings are not passed as a simple pointer, as one might naively guess.  
There's a 2nd hidden integer argument related to the character string length.
Think maybe I got this problem solved now....
*/

data_id *id_new;
i4 nlab,i;
char *fname;
char *datatype;
char *label;

/* 1st, print out the sizes of f77 char arrays (debug) */

/*
printf("len_fnamef77 = %d\n",len_fnamef77);
printf("len_labf77 = %d\n",len_labf77);
printf("len_dtf77 = %d\n",len_dtf77);
printf("len_dataf77 = %d\n",len_dataf77);
*/

/* Convert fortran strings fnamef77,labf77,dtf77 to C-strings fname,label,
   datatype */

fs2cs(&fname,fnamef77,len_fnamef77);
fs2cs(&label,labf77,len_labf77);
fs2cs(&datatype,dtf77,len_dtf77);

/* Debug - print out C-strings fname,label and character variable datatype */

/*
printf("fname = %s\n",fname);
printf("label = %s\n",label);
printf("datatype = %c\n",*datatype);
*/

id_new=malloc(sizeof(data_id));
nlab=strlen(label);
id_new->label=calloc(nlab+1,1);

for (i=0;i<nlab;i++)
{
	*(id_new->label+i)=*(label+i);
}

id_new->nbpw=(i4) *nbpwf77;
id_new->datatype=*datatype;
id_new->order=0; /* value shouldn't matter */
id_new->ndim= (i4) *ndimf77;
id_new->dims=malloc(id_new->ndim * sizeof(pos));

/* I wonder if the f77 integer values of the dimensions shouldn't be
   assumed to be of regular integer type rather than (pos).  
   They could be recast as (pos) inside this fn.  OK, now implemented. */

for (i=0;i<id_new->ndim;i++)
{
	*(id_new->dims+i)= (pos) *(dimsf77+id_new->ndim-1-i); /*flip dims f77*/
}
sdf_replace(fname, (i4) *replace, id_new, dataf77);
free(label);
free(fname);
free(datatype);
free(id_new->label);
free(id_new->dims);
free(id_new);
}

i4 sdf_write(char *fname, data_id *id, void *data)
{

/* Write the data id header and the data into file fname */

/* File structure: (assuming the id-string is "SDF format"):

 * 1.  Bytes 0-10 - file identifier string, including null char terminator 
 * 2.  Bytes 11-18 - 64-bit integer, containing file position of next byte
 *     in the header.  Byte 11 is also strlen(sdfstr)+1 .
 * 3.  Bytes 19-26 - 64-bit integer, containing file position of next byte
       in the data portion of the file (which starts at byte 2019).
       Byte 19 is also strlen(sdfstr)+1+sizeof(pos) .
 * 4.  Bytes 27-30 - 32-bit integer, containing the number of datasets 
       currently written in the file.  Byte 27 is also
       strlen(sdfstr)+1 +2*sizeof(pos) .
 * 5.  Bytes 31-38 - 64-bit integer, containing the size of the file header,
       'hdrsize', counting from byte 19.  Initially this is set to 2000, 
        but can be increased in blocks of 2000 bytes as needed.
 * 6.  Bytes 39-2018 - A series of linefeed terminated strings, containing
       descriptions of the data to be written.  These strings generated by
       the data_to_string function.  Byte 39 is also
       strlen(sdfstr)+1 + 2*sizeof(pos) + sizeof(i4) *sizeof(pos) .
 * 7.  Bytes 2019-eof - the data, with each dataset written squentially,
       in order.  Byte 2019 is also equal to
       strlen(sdfstr)+1+sizeof(pos)+hdrsize .  If hdrsize is larger than
       2000, the data will start correspondingly later in the file
       
       The integers in the header and the data are all written in large-endian
       byte order, if the data isn't simply byte data. */

FILE *fp;
pos hdrpos, datapos, hdrtmp, datatmp, datasize, hdrsize;
pos hdrsztmp, oldhdrsize, startloc;
char *sdfstr="SDF format";
char *filler;
char testid[300];
i4 i,ise,ibe, next_order, next_tmp;
pos hdrinit=2000, safety=100;
int fctest, fotest;
/* fotest is a book-keeping variable that is 1 if file is open, 0 otherwise */
fotest=0;

/* Make certain (pos) is defined as 64bit int, otherwise, exit */

test_sizes(); /* NEW WAY */

/*  OLD WAY
if(sizeof(pos) != 8)
{
	printf("sdf_write: error building sdf_write: sizeof(pos) = %d not 8\n",
            sizeof(pos));
	fflush(stdout);
	exit(1);
}
*/

ibe = is_big_endian ();
ise = 0;
if (ibe == 0) ise = 1;    /* set flag for byteswapping if small endian */

/* Initially, open the file as read-only */

fp=fopen(fname,"rb");

if (fp == NULL) /* does file fname exist? */
{
	/* if file doesnt exist, must create it and initialize it */

/*
        printf("sdf_write:  File %s does not exist; open and initialize\n",
              fname);
*/
        fp=fopen(fname,"wb");
        fotest=1;
        /* write the identifying string 1st thing in the file */
        fwrite(sdfstr,sizeof(char),strlen(sdfstr)+1,fp);
        /* now calculate where initial header info for data will be,
         * after leaving enough space to write in the positions of
         * the header and data parts of the file */
        hdrpos = (pos) ftello(fp)+(pos) 2* sizeof(pos)+(pos) sizeof(i4)
                 +(pos) sizeof(pos);
        hdrtmp=hdrpos;
        /* byteswap the hdr pos if necessary, and write it to file */
	if(ise) byteswap((void *)&hdrtmp,(pos)1,sizeof(pos));
        fwrite(&hdrtmp, sizeof(pos), 1, fp);
        /* now write out 2KB of zeroes, after which the data part will start */
        hdrsize=hdrinit;
        filler=calloc(hdrsize,sizeof(char));
        fwrite(filler,sizeof(char),hdrsize,fp);
        free(filler);
        /* now write the position of the data part to the file */
        datapos = (pos) ftello(fp);
        datatmp=datapos;
        if(ise) byteswap((void *)&datatmp,(pos)1,sizeof(pos));
        fseeko(fp, hdrpos-(pos) sizeof(pos)-(pos)(sizeof(i4))-
                (pos)(sizeof(pos)) ,SEEK_SET); 
        fwrite(&datatmp, sizeof(pos), 1, fp);
        /* now write the order number for next data to be written (0) */
        next_order=(i4) 0;
        next_tmp=next_order;
        /* byte swapping really not necessary since it's 0, but too confusing
         * to not do it */
        if(ise) byteswap((void *)&next_tmp,(pos)1,sizeof(i4));
	fwrite(&next_tmp, sizeof(i4),1,fp);
        hdrsztmp=hdrsize;
        if(ise) byteswap((void *)&hdrsztmp,(pos)1,sizeof(pos));
        fwrite(&hdrsztmp, sizeof(pos),1,fp);
        /* file is now initialized; close it */
        fclose(fp);
	fotest=0;
}
else
{
	fotest=1;
	fgets(testid,strlen(sdfstr)+1,fp);
	/* Now, test to make sure the SDF file identifier string is there */
	if(strncmp(sdfstr,testid,strlen(sdfstr)))
	{
		fclose(fp);
		fotest=0;
		printf("sdf_write: File %s is not an SDF file, exiting\n",
			fname);
		fflush(stdout);
		exit(1);
	}
}

if (fotest) 
{
	fclose(fp); 
	fotest=0;
}

/* re-open file with both read/write permissions */

fp=fopen(fname,"rb+");
fotest=1;

if(fp == NULL)
{
	fotest=0;
	printf("sdf_write:  Can't open file %s\n",fname);
	printf("sdf_write: errno = %d\n",errno);
	fflush(stdout);
	exit(1);
}

/* Now, read the locations for next writes of header info and data */

fseeko(fp, (pos) strlen(sdfstr)+1, SEEK_SET);
fread(&hdrpos,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrpos,(pos)1,sizeof(pos));

fread(&datapos,sizeof(pos),1,fp);
if(ise) byteswap((void *)&datapos,(pos)1,sizeof(pos));

fread(&next_order,sizeof(i4),1,fp);
if(ise) byteswap((void *)&next_order,(pos)1,sizeof(i4));

fread(&hdrsize,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrsize,(pos)1,sizeof(pos));

/* make sure order value written out is actually the order of
 * data written to file */

id->order=next_order;

/* Check to see if header almost full; if so, print warning and increment,
   and record new value at its proper location in the file */

if( hdrpos >= (pos) (hdrsize - safety) )
{
/* now increase size of header area, and move the rest of the file up,
   adjust the current value of datapos accordingly */

 /* For now, comment out this next print statement
 printf("sdf_write: header area almost full, incrementing by hdrinit.\n");  
 fflush(stdout);
 */
 oldhdrsize=hdrsize;
 hdrsize+=hdrinit;
 hdrsztmp=hdrsize;
/*
 fseeko(fp,hdrpos-(pos)(sizeof(pos)),SEEK_SET);
*/
 fseeko(fp,(pos)(strlen(sdfstr)+1+2*sizeof(pos)+sizeof(i4)),SEEK_SET);
 if(ise) byteswap((void *)&hdrsztmp,(pos)1,sizeof(pos));
 fwrite(&hdrsztmp, sizeof(pos),1,fp);
 startloc=(pos) strlen(sdfstr)+1+sizeof(pos)+oldhdrsize;
 move_up(fp,hdrinit,startloc,datapos);
 fflush(fp);
/*
 fn=fileno(fp);
*/
 datapos += hdrinit;
 sdf_file_truncate(fp, (pos) datapos);
/*
 ftruncate(fn, (pos) datapos);
*/
}

/* check that id->datatype is printable, otherwise in big trouble and must
 * exit.  Actually is really only supposed to be f,i, or b acc. to sdf
 * manifesto, but we'll let folks slide here and use anything as long
 * as it's not blank. */

if (!isgraph(id->datatype))
{
	printf("sdf_write: error - id->datatype = 0x%02x not printable\n",
		id->datatype);
	fflush(stdout);
	fclose(fp);
	fotest=0;
	exit(1);
}

/* Now make sure that id->label is 'blank-free': */

rm_blanks(id->label);

/* Now, output the header information for the data to the file header area */

fseeko(fp,hdrpos,SEEK_SET);
data_to_string(id,fp);
hdrpos=ftello(fp);

/* Now, figure out the size of the data file from the info in the id struc: 
 * (Basically, just multiply all the dimensions together) 
 *
 * Note that this was done before the data_size fn was written, so prob.
 * should replace the next 5 lines of code with call to data_size.
*/


datasize=(pos) *(id->dims+0);  /* there's always at least one value */
for (i=1; i<id->ndim; i++)     /* sure hope zero trip loops work */
{
	datasize *=  (pos)( *(id->dims+i) );
}
if(id->datatype != 'c') /* test for complex var */
{
	if(ise) byteswap((void *)data, (pos) datasize, id->nbpw);
}
else
{
	/* kludge for complex var */
	if(ise) byteswap((void *)data, (pos) 2*datasize, id->nbpw);
}

/* Finally, write out the data: */

fseeko(fp,datapos,SEEK_SET);

if (id->datatype != 'c') /* test for complex variable type */
{
	fwrite(data,id->nbpw,datasize,fp); /* is this 64-bit compatible? */
	/* for non-destructive, should btye-swap back! */
	if(ise) byteswap((void *)data, (pos) datasize, id->nbpw);
}
else
{
	fwrite(data,id->nbpw,(pos)2*datasize,fp); /*kludge for complex var */
	/* for non-destructive, should btye-swap back! */
	if(ise) byteswap((void *)data, (pos) 2*datasize, id->nbpw);
}

datapos=(pos) ftello(fp);

/* Finish up: increment next_order, and update hdrpos and datapos in the file */

next_order++;
hdrtmp=hdrpos;
datatmp=datapos;
next_tmp=next_order;
hdrsztmp=hdrsize;

if(ise) byteswap((void *)&hdrtmp,(pos)1,sizeof(pos));
if(ise) byteswap((void *)&datatmp,(pos)1,sizeof(pos));
if(ise) byteswap((void *)&next_tmp,(pos)1,sizeof(i4));
if(ise) byteswap((void *)&hdrsztmp,(pos)1,sizeof(pos));

fseeko(fp, (pos) strlen(sdfstr)+1, SEEK_SET);

fwrite(&hdrtmp,sizeof(pos),1,fp);
fwrite(&datatmp,sizeof(pos),1,fp);
fwrite(&next_tmp,sizeof(i4),1,fp);
fwrite(&hdrsztmp,sizeof(pos),1,fp);

fflush(fp);
if(fotest)
{
	fctest=fclose(fp);
	fotest=0;
}
if(fctest != 0)
{
	printf("sdf_write:  failure to close file %s\n",fname);
	printf("sdf_write: errno = %d\n",errno);
	fotest=1;
}
 
/* We're done */

return 0;
}

i4 sdf_query(char *fname)
{


FILE *fp;
pos datasize;
char *sdfstr="SDF format";
char temp[301], testid[300];
i4 ise,ibe, iorder, ndatasets;
pos hdrsize;
data_id *id;

/* Make certain (pos) is defined as 64bit int, otherwise, exit */

test_sizes(); /* New Way */

/*  OLD WAY
if(sizeof(pos) != 8)
{
	printf("sdf_query: error building sdf_query: sizeof(pos) = %d not 8\n",
            sizeof(pos));
	exit(1);
}
*/

ibe = is_big_endian ();
ise = 0;
if (ibe == 0) ise = 1;    /* set flag for byteswapping if small endian */

/* Initially, open the file as read-only */

fp=fopen(fname,"rb");

if (fp == NULL) /* does file fname exist? */
{
	printf("sdf_query: file %s does not exist, exiting\n",fname);
	fflush(stdout);
	exit(1);
}
else
{
	fgets(testid,strlen(sdfstr)+1,fp);
	/* Now, test to make sure the SDF file identifier string is there */
	if(strncmp(sdfstr,testid,strlen(sdfstr)))
	{
		fclose(fp);
		printf("sdf_query: File %s is not an SDF file, exiting\n",
			fname);
		fflush(stdout);
		exit(1);
	}
}
 
/* read the number of datasets in the file, and the header size: */

fseeko(fp, strlen(sdfstr)+1 + 2*sizeof(pos) , SEEK_SET);
fread(&ndatasets,sizeof(i4),1,fp);
if(ise) byteswap((void *)&ndatasets,(pos)1,sizeof(i4));
fread(&hdrsize,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrsize,(pos)1,sizeof(pos));
/*
printf("ndatasets = %d\n",ndatasets);
*/

/* read the string data of all the datasets, print it out */

for (iorder=0; iorder<ndatasets; iorder++)
{
/*
        fgets(temp,300,fp);
        printf("iorder = %d , header string = %s",iorder,temp);
*/
        id=(data_id *)string_to_data(fp);
/*
	printf("sdf_query debug: id->label =%s",id->label);
*/
	datasize=data_size(id);
	output_int64(temp,datasize);
        printf("sdf_query: dataset %d, size %selements, descriptor: ",
               iorder,temp);
	fflush(stdout);
/* OLD WAY
	#if defined (__SVR4) && defined (__sun) 
        printf("sdf_query: dataset no. %d, datasize = %llu , descriptor = ",
             iorder, datasize);
        #elif defined WIN32
        printf("sdf_query: dataset no. %d, datasize = %I64u , descriptor = ",
             iorder, datasize);
        #else
        printf("sdf_query: dataset no. %d, datasize = %llu , descriptor = ",
             iorder, datasize);
        #endif
*/
        data_to_string(id,stdout);
	fflush(stdout);
	free(id->label);
        free(id->dims);
        free(id);
}

fclose(fp);
return ndatasets;
}

void sdf_housekeeping(char *fname, i4 *norder, pos *hdrposout, 
		pos *dataposout, pos *hdrsizeout)
{

/* return the values of norder, hdrpos, datapos, and hdrsize to the
 * the caller.  Note they are passed by reference, not value.
*/

FILE *fp;
pos hdrpos, datapos;
char *sdfstr="SDF format";
char testid[300];
i4 ise,ibe, ndatasets;
pos hdrsize;

/* Make certain (pos) is defined as 64bit int, otherwise, exit */

test_sizes(); 

ibe = is_big_endian ();
ise = 0;
if (ibe == 0) ise = 1;    /* set flag for byteswapping if small endian */

/* Initially, open the file as read-only */

fp=fopen(fname,"rb");

if (fp == NULL) /* does file fname exist? */
{
	printf("sdf_housekeeping: file %s does not exist, exiting\n",fname);
	fflush(stdout);
	exit(1);
}
else
{
	fgets(testid,strlen(sdfstr)+1,fp);
	/* Now, test to make sure the SDF file identifier string is there */
	if(strncmp(sdfstr,testid,strlen(sdfstr)))
	{
		fclose(fp);
		printf("sdf_housekeeping: File %s is not an SDF file, exiting\n"
			, fname);
		fflush(stdout);
		exit(1);
	}
}
/* read hdrpos,datapos,ndatasets,hdrsize: */
fseeko(fp, strlen(sdfstr)+1, SEEK_SET);

fread(&hdrpos,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrpos,(pos)1,sizeof(pos));

fread(&datapos,sizeof(pos),1,fp);
if(ise) byteswap((void *)&datapos,(pos)1,sizeof(pos));

fread(&ndatasets,sizeof(i4),1,fp);
if(ise) byteswap((void *)&ndatasets,(pos)1,sizeof(i4));

fread(&hdrsize,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrsize,(pos)1,sizeof(pos));

/*
 printf("sdf_housekeeping: existing datasets = %d\n",ndatasets);
 fflush(stdout);
*/

*norder=ndatasets;

/* debug */
/*
output_int64(temp,hdrsize);
printf("sdf_housekeeping: hdrsize = %s\n",temp);
fflush(stdout);
*/

/*
 * Now set output variables to the ones read in from file:
*/
*norder=ndatasets;
*hdrposout = hdrpos;
*dataposout = datapos;
*hdrsizeout = hdrsize;

fclose(fp);
return;
}

void sdf_labmatch_f77__(char *fname_f77, intf * ndat_f77,
     char *labtest_f77, intf *match, intf len_fnamef77, intf len_labtestf77)
{
sdf_labmatch_f77(fname_f77, ndat_f77, labtest_f77, match, len_fnamef77,
    len_labtestf77);
return;
}
void sdf_labmatch_f77_(char *fname_f77, intf * ndat_f77,
     char *labtest_f77, intf *match, intf len_fnamef77, intf len_labtestf77)
{
sdf_labmatch_f77(fname_f77, ndat_f77, labtest_f77, match, len_fnamef77,
    len_labtestf77);
return;
}

void sdf_labmatch_f77(char *fname_f77, intf * ndat_f77,
     char *labtest_f77, intf *match, intf len_fnamef77, intf len_labtestf77)
{
i4 i;
i4 ndat_temp;
i4 *matchind;
char *fname, *labtest;

/* DEBUG - print out lengths of fortran stings (hidden args) */
/*
printf("len_fnamef77 = %d\n",len_fnamef77);
printf("len_labf77 = %d\n",len_labf77);
printf("len_dtf77 = %d\n",len_dtf77);
*/

/* Convert fnamef77 to C-string fname */
fs2cs(&fname,fname_f77,len_fnamef77);
fs2cs(&labtest,labtest_f77,len_labtestf77);
/* DEBUG - print out C-string version of filename fname */
/*
printf("fname = %s\n",fname);
*/
matchind = sdf_labmatch(fname,&ndat_temp,labtest);
*ndat_f77= (intf) ndat_temp;
for (i=0;i<ndat_temp;i++)
{
	*(match + i) = (intf) *(matchind+i);
}

free(fname);
free(labtest);
free(matchind);
return;
}

i4 * sdf_labmatch(char *fname, i4 *ndat, char * labtest)
{
/*
Returns a pointer to an array of dataset indices for which 
id->label matches labtest.  
ndat is determined independently, and its value set in *ndat.  Note that in
the caller, ndat must be allocated before call to sdf_labmatch, or
alternatively its address can be passed to sdf_labmatch.

This function allocates and returns the space for the integer array of
matched indices, and is always of size ndat.  It will be filled with 
values of -1 after the end of matched indices reached.  Be sure to free
the returned index array in caller when finished with it.

*/

FILE *fp;
char *sdfstr="SDF format";
char testid[300];
i4 ise,ibe, iorder, ndatasets;
i4 match, nclabtest, nclab, ncidlab, *matchind;
pos hdrsize;
data_id *id;

/* Make certain (pos) is defined as 64bit int, otherwise, exit */

test_sizes(); /* New way */

/*  OLD WAY
if(sizeof(pos) != 8)
{
	printf("sdf_sizes: error building sdf_sizes: sizeof(pos) = %d not 8\n",
            sizeof(pos));
	exit(1);
}
*/

ibe = is_big_endian ();
ise = 0;
if (ibe == 0) ise = 1;    /* set flag for byteswapping if small endian */

/* Initially, open the file as read-only */

fp=fopen(fname,"rb");

if (fp == NULL) /* does file fname exist? */
{
	printf("sdf_labmatch: file %s does not exist, exiting\n",fname);
	fflush(stdout);
	exit(1);
}
else
{
	fgets(testid,strlen(sdfstr)+1,fp);
	/* Now, test to make sure the SDF file identifier string is there */
	if(strncmp(sdfstr,testid,strlen(sdfstr)))
	{
		fclose(fp);
		printf("sdf_sizes: File %s is not an SDF file, exiting\n",
			fname);
		fflush(stdout);
		exit(1);
	}
}
 
/* read the number of datasets in the file: */

fseeko(fp, strlen(sdfstr)+1 + 2*sizeof(pos) , SEEK_SET);
fread(&ndatasets,sizeof(i4),1,fp);
if(ise) byteswap((void *)&ndatasets,(pos)1,sizeof(i4));
fread(&hdrsize,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrsize,(pos)1,sizeof(pos));

/* set *ndat to ndatasets */
*ndat = ndatasets;

/*
 DEBUG - printf("sdf_size: *ndat = %d\n",*ndat);
*/

/* Now, allocate the memory for *manystrings and matchind */

matchind=malloc(ndatasets*sizeof(i4));
for (iorder=0; iorder<ndatasets; iorder++)
{
	*(matchind+iorder)=(i4) (-1);
}

match=0;
nclabtest=strlen(labtest);
for (iorder=0; iorder<ndatasets; iorder++)
{
        id=(data_id *)string_to_data(fp);
/*
	  printf("sdf_query debug: id->label =%s",id->label);
*/
        ncidlab=strlen(id->label);
        nclab = (nclabtest > ncidlab) ? nclabtest : ncidlab;
        if(!strncmp(id->label,labtest,nclab))
	{
		*(matchind+match)=iorder;
		match++;
	}
	
	free(id->label);
        free(id->dims);
        free(id);
}

fclose(fp);
return matchind;
}
i4 sdf_sizes(char *fname, i4 *ndat, char ** datatypes, pos ** datasizes,
	i4 **nbpw)
{
/*
Returns *ndat, the number of datasets, the char array **datatypes,
the pos array **datasizes, and the i4 array ** nbpw, reflecting the 
*ndat values of datatypes, datasizes, and nbpw values from each dataset
in the sdf file.  Note **datatypes is a true character array, not a string!
After calling sdf_sizes, calling routine responsible for freeing malloc'ed
arrays **datatypes, **datasizes, **nbpw .
Note also that ndat must be allocated BEFORE this function is called.
*/

FILE *fp;
char *sdfstr="SDF format";
char testid[300];
i4 ise,ibe, iorder, ndatasets;
pos hdrsize;
data_id *id;

/* Make certain (pos) is defined as 64bit int, otherwise, exit */

test_sizes(); /* New way */

/*  OLD WAY
if(sizeof(pos) != 8)
{
	printf("sdf_sizes: error building sdf_sizes: sizeof(pos) = %d not 8\n",
            sizeof(pos));
	exit(1);
}
*/

ibe = is_big_endian ();
ise = 0;
if (ibe == 0) ise = 1;    /* set flag for byteswapping if small endian */

/* Initially, open the file as read-only */

fp=fopen(fname,"rb");

if (fp == NULL) /* does file fname exist? */
{
	printf("sdf_sizes: file %s does not exist, exiting\n",fname);
	fflush(stdout);
	exit(1);
}
else
{
	fgets(testid,strlen(sdfstr)+1,fp);
	/* Now, test to make sure the SDF file identifier string is there */
	if(strncmp(sdfstr,testid,strlen(sdfstr)))
	{
		fclose(fp);
		printf("sdf_sizes: File %s is not an SDF file, exiting\n",
			fname);
		fflush(stdout);
		exit(1);
	}
}
 
/* read the number of datasets in the file: */

fseeko(fp, strlen(sdfstr)+1 + 2*sizeof(pos) , SEEK_SET);
fread(&ndatasets,sizeof(i4),1,fp);
if(ise) byteswap((void *)&ndatasets,(pos)1,sizeof(i4));
fread(&hdrsize,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrsize,(pos)1,sizeof(pos));

/* set *ndat to ndatasets */
*ndat = ndatasets;

/*
 DEBUG - printf("sdf_size: *ndat = %d\n",*ndat);
*/

/* Now, allocate the memory for *datatypes, *datasizes, *nbpw: */

*datatypes=malloc(ndatasets*sizeof(char));
*datasizes=malloc(ndatasets*sizeof(pos));
*nbpw=malloc(ndatasets*sizeof(i4));

for (iorder=0; iorder<ndatasets; iorder++)
{
/*
        fgets(temp,300,fp);
        printf("iorder = %d , header string = %s",iorder,temp);
*/
        id=(data_id *)string_to_data(fp);
/*
	  printf("sdf_query debug: id->label =%s",id->label);
*/
	*(*datasizes+iorder)=data_size(id);
	*(*datatypes+iorder)=id->datatype;
	*(*nbpw+iorder)=id->nbpw;
	
	free(id->label);
        free(id->dims);
        free(id);
}

fclose(fp);
return 0;
}

void sdf_sizes_f77__(char *fnamef77, intf *ndatf77, char dtypesf77[][1],
intf * datasizesf77, intf *nbpwf77, intf len_fnamef77, intf len_dtypesf77)
{
sdf_sizes_f77(fnamef77, ndatf77, dtypesf77,
datasizesf77, nbpwf77, len_fnamef77, len_dtypesf77);
return;
}
void sdf_sizes_f77_(char *fnamef77, intf *ndatf77, char dtypesf77[][1],
intf * datasizesf77, intf *nbpwf77, intf len_fnamef77, intf len_dtypesf77)
{
sdf_sizes_f77(fnamef77, ndatf77, dtypesf77,
datasizesf77, nbpwf77, len_fnamef77, len_dtypesf77);
return;
}

void sdf_sizes_f77(char *fnamef77, intf *ndatf77, char dtypesf77[][1],
intf * datasizesf77, intf *nbpwf77, intf len_fnamef77, intf len_dtypesf77)
{

/*
Fortran 77 call-able version of sdf_sizes - returns size information on
each dataset in an sdf file 
*/
char *dtypes;
pos *dsizes;
i4 *nbpw;
i4 *ndat;
i4 i;
char *fname;
ndat=malloc(sizeof(i4)); /*note ndat must be allocated before sdf_sizes called*/
fs2cs(&fname,fnamef77,len_fnamef77); /*convert fname from fstring to cstring */
sdf_sizes(fname, ndat, &dtypes, &dsizes, &nbpw); /* get dataset sizes */

/*
 DEBUG - printf("sdf_sizes_f77: *ndat = %d\n",*ndat);
*/

for (i=0;i<*ndat;i++)
{
	*(nbpwf77+i)= (intf) *(nbpw+i);
	*(datasizesf77+i)=(intf) *(dsizes+i); /* converted to intf for f77 */
	dtypesf77[i][0]=(char)*(dtypes+i);
}
*ndatf77= (intf) *ndat;

/* free all the temporary junk created in this function */

free(ndat);
free(fname);
free(dtypes);
free(nbpw);
free(dsizes);
return;
}
void sdf_rm_f77__(char *fnamef77, intf len_fnamef77)
{
sdf_rm_f77(fnamef77, len_fnamef77);
return;
}
void sdf_rm_f77_(char *fnamef77, intf len_fnamef77)
{
sdf_rm_f77(fnamef77, len_fnamef77);
return;
}
void sdf_rm_f77(char *fnamef77, intf len_fnamef77)
{
/* Fortran callable function to remove sdf files */
char *fname;
fs2cs(&fname,fnamef77,len_fnamef77);
sdf_rm(fname);
free(fname);
return;
}

i4 sdf_rm(char *fname)
/*
Function to remove the file fname, but only if it's an sdf file.
*/
{
FILE *fp;
char *sdfstr="SDF format";
char testid[300];

fp=fopen(fname,"rb");

if (fp == NULL) /* does file fname exist? */
{
	printf("sdf_rm: file %s cannot be opened for deletion, returning\n",
          fname);
	fflush(stdout);
	return 1;
}

fgets(testid,strlen(sdfstr)+1,fp);
/* Now, test to make sure the SDF file identifier string is there */
if(strncmp(sdfstr,testid,strlen(sdfstr)))
{
	fclose(fp);
	printf("sdf_rm: File %s is not an SDF file, cannot remove\n",
		fname);
	fflush(stdout);
	exit(1);
}
fclose(fp);
remove(fname);
return 0;
}

void sdf_query_f77__(char *fnamef77, intf *ndat, intf len_fname_f77)
{
   sdf_query_f77(fnamef77, ndat, len_fname_f77);
   return;
}
void sdf_query_f77_(char *fnamef77, intf *ndat, intf len_fname_f77)
{
   sdf_query_f77(fnamef77, ndat, len_fname_f77);
   return;
}
void  sdf_query_f77(char *fnamef77, intf * ndat, intf len_fname_f77)
{

/* This is the fortran callable version of sdf_query. *ndat is
the number of datasets in the file. */

char *sdfstr="SDF format";
char *fname;
i4 ise,ibe;

/* Make certain (pos) is defined as 64bit int, otherwise, exit */

test_sizes(); /* New way */

/* OLD WAY
if(sizeof(pos) != 8)
{
   printf("sdf_query_f77: error building, sizeof(pos) = %d not 8\n",
            sizeof(pos));
   exit(1);
}
*/

ibe = is_big_endian ();
ise = 0;
if (ibe == 0) ise = 1;    /* set flag for byteswapping if small endian */

/*
OK, now convert the Fortran string containing the name of the file to the
C-string containing the name of the file, and call sdf_query:
*/

fs2cs(&fname, fnamef77, len_fname_f77);
*ndat=sdf_query(fname);
free(fname);
return;
}

data_id *sdf_read(char *fname, i4 order, void **data)
{

/* read in the dataset corresponding to order, and return idnew */

FILE *fp;
data_id *id;
pos datasize, curpos;
char *sdfstr="SDF format";
char temp[301], testid[300];
i4 ise,ibe, iorder, ndatasets;
pos hdrsize;

/* Make certain (pos) is defined as 64bit int, otherwise, exit */

test_sizes(); /* New way */

/* OLD WAY
if(sizeof(pos) != 8)
{
	printf("sdf_read: error building sdf_read: sizeof(pos) = %d not 8\n",
            sizeof(pos));
	exit(1);
}
*/

ibe = is_big_endian ();
ise = 0;
if (ibe == 0) ise = 1;    /* set flag for byteswapping if small endian */

/* Initially, open the file as read-only */

fp=fopen(fname,"rb");

if (fp == NULL) /* does file fname exist? */
{
	printf("sdf_read: file %s does not exist, exiting\n",fname);
	fflush(stdout);
	exit(1);
}
else
{
	fgets(testid,strlen(sdfstr)+1,fp);
	/* Now, test to make sure the SDF file identifier string is there */
	if(strncmp(sdfstr,testid,strlen(sdfstr)))
	{
		fclose(fp);
		printf("sdf_read: File %s is not an SDF file, exiting\n",
			fname);
		fflush(stdout);
		exit(1);
	}
}
/* read the number of datasets in the file: */

fseeko(fp, strlen(sdfstr)+1 + 2*sizeof(pos) , SEEK_SET);
fread(&ndatasets,sizeof(i4),1,fp);
if(ise) byteswap((void *)&ndatasets,(pos)1,sizeof(i4));
/*
 Next print statement possibly important but nuke for now
*/

/*
printf("sdf_read: existing datasets = %d\n",ndatasets);
fflush(stdout);
*/

/*
printf("sdf_read: desired dataset order = %d\n",order);
*/
if(order > (ndatasets-1))
{
	printf("sdf_read:  Desired dataset order %d > ndatasets-1 = %d\n",
             order, ndatasets-1);
	fflush(stdout);
	exit(1);
}
/* read in header size */
fread(&hdrsize,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrsize,(pos)1,sizeof(pos));
/* debug */
/*
output_int64(temp,hdrsize);
printf("sdf_read: hdrsize = %s\n",temp);
fflush(stdout);
*/

/* read the string data of all the datasets, print it out */
curpos = (pos) strlen(sdfstr)+1+sizeof(pos)+hdrsize;
/*
        printf("DEBUG - sdf_read: curpos = %qu\n",curpos);
*/
for (iorder=0; iorder <= order; iorder++)
{
/*
        fgets(temp,300,fp);
        printf("iorder = %d , header string = %s",iorder,temp);
*/

        id=string_to_data(fp);
	datasize=data_size(id);


        if(iorder < order)
	{
		if(id->datatype != 'c') /* test for complex var */
		{
			curpos += (pos) datasize*(id->nbpw);
		}
		else
		{
			curpos += (pos) datasize*2*(id->nbpw);
			/* kludge for complex variables */
		}
		free(id->label);
		free(id->dims);
		free(id);
	}
	else
	{
        /* NEW WAY */
	output_int64(temp,datasize);
	/* next print statements possibly important, but nuke for now */
        /*
        printf("sdf_read: dataset %d, size %selements, descriptor: ",
               iorder,temp);
	fflush(stdout);
        */

        /* OLD WAY 
	#if defined (__SVR4) && defined (__sun) 
        printf("sdf_read: dataset no. %d, datasize = %llu , descriptor = ",
             iorder, datasize);
        #elif defined WIN32
        printf("sdf_read: dataset no. %d, datasize = %I64u , descriptor = ",
             iorder, datasize);
        #else
        printf("sdf_read: dataset no. %d, datasize = %llu , descriptor = ",
             iorder, datasize);
        #endif
        */
        /*
        data_to_string(id,stdout);
	printf("\n");
	fflush(stdout);
        */
	}
}

fseeko(fp,curpos,SEEK_SET);                /* reposition file pointer  */

if (id->datatype != 'c') /* test for complex var */
{
	*data=malloc(datasize* (pos) (id->nbpw));         /* allocate data */
	fread(*data, id->nbpw , datasize,fp);      /* read the data */
	if(ise) byteswap(*data,datasize,id->nbpw); /* byteswap if neeeded */
}
else
{
	/* kludge for reading complex variable data */
	*data=malloc(datasize*(pos)2*(id->nbpw)); 
	fread(*data, id->nbpw , 2*datasize,fp);      /* read the data */
	if(ise) byteswap(*data,2*datasize,id->nbpw); /* byteswap if neeeded */
}

fclose(fp);
return id;

}

void sdf_details_f77__(char *fnamef77, intf *iorder, char *labf77, char *dtf77, 
     intf *nbpwf77, intf *ndimf77, intf *dimsf77,
     intf len_fnamef77, intf len_labf77, intf len_dtf77)
{
sdf_details_f77(fnamef77, iorder, labf77, dtf77, 
     nbpwf77, ndimf77, dimsf77,
     len_fnamef77, len_labf77, len_dtf77);
return;
}
void sdf_details_f77_(char *fnamef77, intf *iorder, char *labf77, char *dtf77, 
     intf *nbpwf77, intf *ndimf77, intf *dimsf77,
     intf len_fnamef77, intf len_labf77, intf len_dtf77)
{
sdf_details_f77(fnamef77, iorder, labf77, dtf77, 
     nbpwf77, ndimf77, dimsf77,
     len_fnamef77, len_labf77, len_dtf77);
return;
}

void sdf_details_f77(char *fnamef77, intf *iorder, char *labf77, char *dtf77, 
     intf *nbpwf77, intf *ndimf77, intf *dimsf77,
     intf len_fnamef77, intf len_labf77, intf len_dtf77)

{
data_id *id;
i4 nlab,i;
char *fname;
pos datasize;

/* DEBUG - print out lengths of fortran stings (hidden args) */
/*
printf("len_fnamef77 = %d\n",len_fnamef77);
printf("len_labf77 = %d\n",len_labf77);
printf("len_dtf77 = %d\n",len_dtf77);
*/

/* Convert fnamef77 to C-string fname */
fs2cs(&fname,fnamef77,len_fnamef77);
/* DEBUG - print out C-string version of filename fname */
/*
printf("fname = %s\n",fname);
*/

id=(data_id *) sdf_details(fname, (i4) *iorder);
/*
 id=(data_id *) sdf_read(fname, *iorder, (void **)&dataf77);
*/
nlab=strlen(id->label);

/* now convert id->label into fortran string labf77 */
cs2fs(id->label,labf77,len_labf77);

*dtf77 = id->datatype;
*nbpwf77 = (intf) id->nbpw;
*ndimf77 = (intf) id->ndim;
for (i=0; i < id->ndim; i++)
{
	/* flip order for f77 */
	*(dimsf77+i) = (intf) *(id->dims+id->ndim-1 -i); 
}
datasize=data_size(id);

/* rather than go through the next intermediate stage, could call sdf_read with
 * dataf77 directly, maybe?  ...tried that, didn't work, but prob. needs
 * further investigation ... */

free(fname);
free(id->label);
free(id->dims);
free(id);
}
data_id *sdf_details(char *fname, i4 order)
{

/* return the id structure corresponding to dataset iorder */

FILE *fp;
data_id *id;
pos datasize;
char *sdfstr="SDF format";
char temp[301], testid[300];
i4 ise,ibe, iorder, ndatasets;
pos hdrsize;

/* Make certain (pos) is defined as 64bit int, otherwise, exit */

test_sizes(); 

ibe = is_big_endian ();
ise = 0;
if (ibe == 0) ise = 1;    /* set flag for byteswapping if small endian */

/* Initially, open the file as read-only */

fp=fopen(fname,"rb");

if (fp == NULL) /* does file fname exist? */
{
	printf("sdf_details: file %s does not exist, exiting\n",fname);
	fflush(stdout);
	exit(1);
}
else
{
	fgets(testid,strlen(sdfstr)+1,fp);
	/* Now, test to make sure the SDF file identifier string is there */
	if(strncmp(sdfstr,testid,strlen(sdfstr)))
	{
		fclose(fp);
		printf("sdf_details: File %s is not an SDF file, exiting\n",
			fname);
		fflush(stdout);
		exit(1);
	}
}
/* read the number of datasets in the file: */

fseeko(fp, strlen(sdfstr)+1 + 2*sizeof(pos) , SEEK_SET);
fread(&ndatasets,sizeof(i4),1,fp);
if(ise) byteswap((void *)&ndatasets,(pos)1,sizeof(i4));
printf("sdf_details: existing datasets = %d\n",ndatasets);
fflush(stdout);

/*
printf("sdf_details: desired dataset order = %d\n",order);
*/
if(order > (ndatasets-1))
{
	printf("sdf_read:  Desired dataset order %d > ndatasets-1 = %d\n",
             order, ndatasets-1);
	fflush(stdout);
	exit(1);
}
/* read in header size */
fread(&hdrsize,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrsize,(pos)1,sizeof(pos));
/* debug */
/*
output_int64(temp,hdrsize);
printf("sdf_details: hdrsize = %s\n",temp);
fflush(stdout);
*/

/* read the string data of all the datasets, print it out */

for (iorder=0; iorder <= order; iorder++)
{
/*
        fgets(temp,300,fp);
        printf("iorder = %d , header string = %s",iorder,temp);
*/

        id=string_to_data(fp);
	datasize=data_size(id);


        if(iorder < order)

	{
		free(id->label);
		free(id->dims);
		free(id);
	}

	else

	{
	output_int64(temp,datasize);
        printf("sdf_details: dataset %d, size %selements, descriptor: ",
               iorder,temp);
	fflush(stdout);
        data_to_string(id,stdout);
	printf("\n");
	fflush(stdout);
	}
}

fclose(fp);
return id;
}

data_id * string_to_data( FILE *fp )
{

/* Read in header string, extract information, return it as a data_id struc */

	char temp[301];
	int ipos,i,ii;
	char *p, *last;
	char * delim = " ";
        data_id *idnew;

        idnew=malloc(sizeof(data_id)); /* create space for structure idnew */
        fgets(temp, 300, fp);          /* read in header string */

/*
	printf("DEBUG - string_to_data: temp = %s\n",temp);
*/

	p=(char *)strtok_r(temp,delim,&last);  /* get 1st token from string */

        /* loop over remaining tokens in the string and put them into idnew */

        i=0;
        ipos=0;
	idnew->order=atol(p);
	while ((p = (char *)strtok_r(NULL,delim,&last)) != NULL)
	{
                i++;
/*
		printf(" string_to_data DEBUG - i = %u, token = %s\n",i,p);
		printf("strlen(p) = %d\n",strlen(p));
*/
		if(i == 1) 
		{
			idnew->label=calloc(strlen(p)+1,1);
/*
			strncpy(idnew->label,p,strlen(p));
*/
			for (ii=0;ii<=strlen(p);ii++)
			{
                		*(idnew->label+ii)=*(p+ii);
			}
		}
                if(i == 2) idnew->datatype=*p;
                if(i == 3) idnew->nbpw=atol(p);
                if(i == 4) idnew->ndim=atol(p);

		/* now create space for the individual dimension values */

		if (i == 5) idnew->dims=malloc(idnew->ndim*sizeof(pos));

                if (i >= 5)
		{
			if(ipos < idnew->ndim)
			{
				/* read in the array dimensions */
				*(idnew->dims+ipos)= (pos) atopos(p);	
				ipos++;
			}
		}
	}
	return idnew;
}
int data_to_buff(data_id *id, char *buff, pos *offset)
{

/* take information from the id structure, and write it as string
to the buffer buff; return the new offset.  This function very similar
to data_to_string function, which writes the string directly to file */

        char *delim = " ";
        char *na = "NA";
        /* char *p, *label; */   /* not used anymore */
	pos ipos;
        i4 i,blank_label,ic,lablen;
        ipos=0;
        sprintf((buff+ipos),"%u ",id->order);
        ipos=(pos) strlen(buff);

	/* Here, need to check that label is not a blank or null string; 
	if so, set to "NA". */

        blank_label=1;
	lablen=strlen(id->label);
        for (ic=0; ic < lablen; ic++)
	{
		if( isprint(*(id->label+ic)) && ! isspace(*(id->label+ic)) )
                     blank_label=0;
	}
/*
	if(blank_label) id->label="NA";
*/
	if(blank_label) strncpy(id->label,na,2);

	sprintf((buff+ipos),"%s ",id->label);

/*
	printf("inside data_to_buff, id->label = %s\n",id->label);
*/

	ipos=(pos) strlen(buff);
        sprintf((buff+ipos),"%c ",id->datatype);
        ipos=(pos) strlen(buff);
        sprintf((buff+ipos),"%u ",id->nbpw);
        ipos=(pos) strlen(buff);
        sprintf((buff+ipos),"%u ",id->ndim);
        for (i=0;i<id->ndim;i++)
	{
		ipos=(pos) strlen(buff);
		output_int64((buff+ipos),*(id->dims+i)); /* new way */
	}

        ipos=(pos) strlen(buff);
        *(buff+ipos)='\n'; /* write linefeed to buffer */
        ipos++;
/*
        printf("data_to_buff: ipos = %d\n",ipos);
*/
        *offset=ipos;
        return 0;
}
int data_to_string(data_id *id, FILE *fp)
{

/* take information from the id structure, and write it to header as string */

	char temp[300];
        char *delim = " ";
        char *na = "NA";
        /* char *p, *label; */  /* not used anymore */
        i4 ipos,i,blank_label,ic,lablen;
        ipos=0;
        sprintf((temp+ipos),"%u ",id->order);
        ipos=strlen(temp);

	/* Here, need to check that label is not a blank or null string; 
	if so, set to "NA". */

        blank_label=1;
	lablen=strlen(id->label);
        for (ic=0; ic < lablen; ic++)
	{
		if( isprint(*(id->label+ic)) && ! isspace(*(id->label+ic)) )
                     blank_label=0;
	}
/*
	if(blank_label) id->label="NA";
*/
	if(blank_label) strncpy(id->label,na,2);

	sprintf((temp+ipos),"%s ",id->label);

/*
	printf("inside data_to_string, id->label = %s\n",id->label);
*/

	ipos=strlen(temp);
        sprintf((temp+ipos),"%c ",id->datatype);
        ipos=strlen(temp);
        sprintf((temp+ipos),"%u ",id->nbpw);
        ipos=strlen(temp);
        sprintf((temp+ipos),"%u ",id->ndim);
        for (i=0;i<id->ndim;i++)
	{
		ipos=strlen(temp);
		output_int64((temp+ipos),*(id->dims+i)); /* new way */

		/*  OLD WAY
		#if defined (__SVR4) && defined (__sun) 
		   sprintf((temp+ipos),"%llu ",*(id->dims + i));
		#elif defined WIN32
	      	   sprintf((temp+ipos),"%I64u ",*(id->dims + i));
		#else
		   sprintf((temp+ipos),"%llu ",*(id->dims + i));
		#endif
		*/
	}
/*
        fprintf(fp, "%s", temp); 
*/
        /* write out the string to the file */
        /* add a line feed */
/*
        fputc((unsigned char)'\n',fp); 
*/

        fprintf(fp, "%s\n", temp); /* write string and linefeed to file */
        return 0;
}

int output_int64(char *buff, pos ill)

/* Function to print a 64-bit integer into the string buff,
   since there seems to be no platform independent format stmt for 64 bit 
   integers.  Adds a trailing blank to the string. User must make sure
   buff is dimensioned or allocated big enough in calling program.  */

{
char work[100];
int neg=0;
int rem,i;
int counter=0;

if(ill < 0)
{
	ill=-ill;
	neg=1;
}

while(ill > 0)
{
	/* digit by digit, get the chars for each decimal place into work */
	rem= (int) (ill % (pos)10);
	ill /= (pos)10;
	sprintf(&work[counter],"%d",rem);
	counter++;
}

if(neg)
{
	work[counter]='-'; /* add - sign for neg. integers */
	counter++;
}

if(counter == 0)
{
	/* make sure that if integer = 0, that at least a 0 gets printed */
	work[counter]='0';
	counter++;
}

counter++; /* add one extra character for the trailing blank */

for (i=0;i<counter-1;i++)
{
	*(buff+counter-i-2)=work[i]; /* digits and minus sign if it exists */
}
        *(buff+counter-1)=' '; /* trailing blank */
	*(buff+counter)='\0'; /* null terminator */

return 0;
}

pos atopos(char *s)
/* Try to make a platform-independent version of atoll, which doesn't work
on Loraine's MAC laptop running OS 10.3.9 .  Returns the value of the
64-bit integer corresponding to the ascii string s. */
{
i4 nchar;
i4 i;
pos value=0;
i4 neg=0;
nchar=strlen(s);
/* should put in an upper limit on nchar */
if(nchar > 19)
{
	printf("atopos: string too long to convert to pos, > 19 chars\n");
	exit(1);
}
for(i=0;i<nchar;i++)
{
	if((s[i] >= '0') && (s[i] <= '9')) 
	{
		value=10*value+(s[i]-'0');
		/* build up the integer value digit by digit */
	}
	if(s[i] == '-') neg=1;
}
if(neg) value=-value;
return value;
}

pos data_size(data_id * id)

/* compute the size of the dataset (in elements) from the id structure */

{
   pos datasize;
   i4 i;
   datasize=(pos) *(id->dims+0);  /* there's always at least one value */
   for (i=1; i<id->ndim; i++)     /* sure hope zero trip loops work */
   {
	datasize *=  (pos)( *(id->dims+i) );
   }
   return datasize;
}

i4 cs2fs(char *cs, char *fs, intf len_fs)
{
/*
  Convert the C string *cs to the Fortran string *fs.  Here it is assumed
  the C string *cs already exists, and the fortran string *fs has already
  been allocated by the fortran calling program, and has length len_fs.
*/
i4 i,csl;
if (len_fs <= 0)
{
	printf("cs2fs:  Bad value of len_fs = %d\n",len_fs);
	fflush(stdout);
	exit(1);
}
csl=strlen(cs);
if ( (csl < 0) || (csl > len_fs))
{
	printf("cs2fs:  Bad value of csl = %d\n",csl);
	fflush(stdout);
	exit(1);
}
for (i=0;i<csl;i++)
{
	*(fs+i)=*(cs+i);
}
for (i=csl;i<len_fs;i++)
{
	*(fs+i)=' ';
}
return 0;
}

i4 fs2cs(char **cs, char *fs, intf len_fs)
{
/*
   Convert the fortran string *fs to the c-style string *cs.  Fortran strings
   are right-filled with blanks; C strings are null terminated.
   Start from end of fortran string, find first non-blank character, then
   copy from the start to that character into the c-style string, add null
   character.  The parameter len_fs has been acquired from the fortran
   call to a C function already.  This malloc's *cs as a new string, so
   must free cs in the calling function when done.
*/
i4 i,csl;
if (len_fs <= 0)
{
	printf("fs2cs:  Bad value of len_fs = %d\n",len_fs);
	fflush(stdout);
	exit(1);
}
csl=len_fs;
for (i=len_fs-1; i >= 0; i--) /* this should be a while loop - cleaner */
{
	if(*(fs+i) != ' ') 
	{
		csl=i+1;
		i=0;
	}
}
/*
 printf("string fs has nontrivial length of %d\n",csl);
*/
*cs=malloc(csl+1);
for (i=0;i<csl;i++)
{
	*(*cs+i)=*(fs+i);
}
*(*cs+csl)='\0';
return 0;
}

void rm_blanks(char *s)
{
/* remove any blanks from a character string.  This is needed for a blank-free
 * string in id->label in sdf_write */
char *work;
i4 i,ls,inb;
ls=strlen(s);
work=malloc(ls+1);
inb=0;
for (i=0;i<ls;i++)
{
	if(*(s+i) != ' ')
	{
		*(work+inb)=*(s+i);
		inb++;
	}
}
*(work+inb)='\0';
for (i=0; i <= inb; i++)
{
	*(s+i)=*(work+i);
}
free(work);
return;
}

i4 move_up(FILE *fp, pos shiftup, pos startloc, pos endloc)
/*

This function shifts the contents of a file upward in the file
by an amount shiftup bytes, starting at position startloc, and ending
at endloc (counted before the shift is done).  endloc is assumed to be the
position at the end of the file before this process starts.
The file (FILE *fp) is assumed to have already been opened.

*/
{
pos bigbuf = 1048576; /* a 1 MB max buffer size, could be adjusted */
pos bufsize;
pos sz_shift;
pos nbuf, lower_loc;
pos fag_end;
char *buf;
pos i;
sz_shift=(endloc-startloc);
nbuf=sz_shift/bigbuf;
fag_end=sz_shift % bigbuf;
bufsize = (nbuf > 0) ? bigbuf : fag_end;
buf = malloc(bufsize);
lower_loc=endloc;
for (i=0;i <= nbuf;i++)
{
	bufsize = (nbuf > i) ? bigbuf : fag_end;
        lower_loc -= bufsize;
	fseeko(fp,lower_loc,SEEK_SET);
/* debug
        printf("lower_loc = %llu\n",lower_loc);
*/
	fread(buf,bufsize,1,fp);
        fseeko(fp,lower_loc+shiftup,SEEK_SET);
        fwrite(buf,bufsize,1,fp);
/* debug 
        location=ftello(fp);
        printf("location = %llu\n",location);
        fflush(stdout);
*/
}

free(buf);
return 0;
}

i4 move_dn(FILE *fp, pos shiftdn, pos startloc, pos endloc)
/*

This function shifts the contents of a file downward in the file
by an amount shiftdn bytes, starting at position startloc, and ending
at endloc (counted before the shift is done).  endloc is assumed to be the
position at the end of the file before this process starts.
The file (FILE *fp) is assumed to have already been opened, and remains open.

*/
{
pos bigbuf = 1048576; /* a 1 MB max buffer size, could be adjusted */

pos bufsize;
pos sz_shift;
pos nbuf, lower_loc;
pos fag_end;
char *buf;
pos i;

/* debug 
printf("shiftdn = %llu, startloc = %llu, endloc = %llu\n",shiftdn,startloc,
          endloc);
*/

sz_shift=(endloc-startloc-shiftdn); /* note this is diff. than in move_up! */
nbuf=sz_shift/bigbuf;
fag_end=sz_shift % bigbuf;
bufsize = (nbuf > 0) ? bigbuf : fag_end;

if(bufsize > 0)
{

   /* debug
   printf("move_dn: bufsize = %llu\n",bufsize);
   */

   buf = malloc(bufsize);
   lower_loc=startloc;
   for (i=0;i <= nbuf;i++)
   {
	bufsize = (nbuf > i) ? bigbuf : fag_end;

   /* debug
        printf("i = %llu\n",i);
        printf("lower_loc = %llu\n",lower_loc);
        printf("bufsize = %llu\n", bufsize);
        printf("shiftdn = %llu\n", shiftdn);
        fflush(stdout);
  */

	fseeko(fp,lower_loc+shiftdn,SEEK_SET);
	fread(buf,bufsize,1,fp);
        fseeko(fp,lower_loc,SEEK_SET);
        fwrite(buf,bufsize,1,fp);
        lower_loc+=bufsize;
   }
   free(buf);

}

/*
If after shifting down, you want to close the file and have it be
shorter, do the following after the call to shift_dn:

fflush(fp);
fn=fileno(fp);
ftruncate(fn, (pos)endloc - (pos)shiftdn )

*/
return 0;
}

i4 sdf_delete(char *fname, i4 idelete)
{

/* Delete dataset no. idelete from the file fname */

FILE *fp;
data_id *id;
pos hdrpos, datapos, hdrtmp, datatmp, datasize, curpos, off, offset;
pos size_del, byte_del, filepos;
char *sdfstr="SDF format";
char *buff;
char testid[300];
i4 ise,ibe, next_order, next_tmp, bufsize, iorder;
pos hdrinit=2000, safety=100;
pos hdrsize, hdrsztmp;
int fotest;

/* Make certain (pos) is defined as 64bit int, otherwise, exit */

test_sizes();
fotest=0;


ibe = is_big_endian ();
ise = 0;
if (ibe == 0) ise = 1;    /* set flag for byteswapping if small endian */

/* Initially, open the file as read-only */

fp=fopen(fname,"rb");

if (fp == NULL) /* does file fname exist? */
{
	/* if not, print msg and exit */
        printf("sdf_delete: file %s does not exist; exiting\n",fname);
	fflush(stdout);
        exit(1);
}
else
{
	fotest=1;
	fgets(testid,strlen(sdfstr)+1,fp);
	/* Now, test to make sure the SDF file identifier string is there */
	/* If not, print msg and exit */
	if(strncmp(sdfstr,testid,strlen(sdfstr)))
	{
		fclose(fp);
		printf("sdf_delete: File %s is not an SDF file, exiting\n",
			fname);
		fflush(stdout);
                fotest=0;
		exit(1);
	}
}

if (fotest)
{
	/*close file if it is open */
        fclose(fp);
        fotest=0;
}

/* re-open file with both read/write permissions */

fp=fopen(fname,"rb+");

/* Now, read the locations for next writes of header info and data */

fseeko(fp, (pos) strlen(sdfstr)+1, SEEK_SET);
fread(&hdrpos,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrpos,(pos)1,sizeof(pos));

fread(&datapos,sizeof(pos),1,fp);
if(ise) byteswap((void *)&datapos,(pos)1,sizeof(pos));

fread(&next_order,sizeof(i4),1,fp);
if(ise) byteswap((void *)&next_order,(pos)1,sizeof(i4));

fread(&hdrsize,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrsize,(pos)1,sizeof(pos));

if(idelete >= next_order)
{
	printf("sdf_delete: can't delete non-existent dataset no. %d in %s\n",
            idelete,fname);
	fflush(stdout);
        fclose(fp);
        exit(1);
}
if(next_order == 1)
{
	fclose(fp);
        fotest=0;
	sdf_rm(fname);
	return 0;
}

bufsize=hdrsize-sizeof(pos)-sizeof(i4); 
buff=(char *)calloc(bufsize*sizeof(char),1);

off=(pos)0;
curpos = (pos) strlen(sdfstr)+1+sizeof(pos)+hdrsize;
for (iorder=0; iorder < next_order; iorder++)
{
	id=(data_id *) string_to_data(fp);
        datasize=data_size(id);
	filepos=ftello(fp);
	if(iorder == idelete)
	{
		size_del = data_size(id);
                byte_del = size_del * id->nbpw;
		if(id->datatype == 'c')
		{
			byte_del *= (pos)2;
		}
		/* Remove the offending dataset, shorten file, reposition */

                /* debug
                printf("byte_del = %d, curpos = %d, datapos = %d\n",
                    (i4)byte_del,(i4)curpos,(i4)datapos);
                */

		move_dn(fp, byte_del, curpos, datapos);
		fflush(fp);
/*
		fn=fileno(fp);
		ftruncate(fn, (pos) datapos - (pos) byte_del );
*/
                sdf_file_truncate(fp, (pos) datapos - (pos) byte_del );
		fseeko(fp,filepos,SEEK_SET);

	}
	if(iorder >= idelete) /* right limit for this? */
	{
		id->order--;
	}
	if(iorder != idelete)
	{
		data_to_buff(id,buff+off,&offset);
		off+=offset;
	}
	if(id->datatype != 'c') /* test for complex var */
        {
                curpos += (pos) datasize*(id->nbpw);
        }
        else
        {
                curpos += (pos) datasize*2*(id->nbpw);
                /* kludge for complex variables */
        }
	free(id->label);
	free(id->dims);
	free(id);
}
hdrpos = strlen(sdfstr)+1 + 2*sizeof(pos) + sizeof(i4) + sizeof(pos);
fseeko(fp, hdrpos, SEEK_SET);
/* Gonna delete the \0 here */
/* fprintf(fp,"%s%c",buff,'\0'); */
fprintf(fp,"%s",buff);
hdrpos= (pos) ftello(fp);
datapos -= byte_del;
next_order--;

hdrtmp=hdrpos;
datatmp=datapos;
next_tmp=next_order;
hdrsztmp=hdrsize;

if(ise) byteswap((void *)&hdrtmp,(pos)1,sizeof(pos));
if(ise) byteswap((void *)&datatmp,(pos)1,sizeof(pos));
if(ise) byteswap((void *)&next_tmp,(pos)1,sizeof(i4));
if(ise) byteswap((void *)&hdrsztmp,(pos)1,sizeof(pos));

fseeko(fp, (pos) strlen(sdfstr)+1, SEEK_SET);

fwrite(&hdrtmp,sizeof(pos),1,fp);
fwrite(&datatmp,sizeof(pos),1,fp);
fwrite(&next_tmp,sizeof(i4),1,fp);
fwrite(&hdrsztmp,sizeof(pos),1,fp);

fclose(fp);
free(buff);

/* We're done */

return 0;
}

i4 sdf_insert(char *fname, i4 insert, data_id * id_new, void * data_new)
{

/* Insert new dataset after no. "insert" from the file fname */

FILE *fp;
data_id *id;
pos hdrpos, datapos, hdrtmp, datatmp, datasize, curpos, off, offset;
pos size_ins, byte_ins, filepos;
char *sdfstr="SDF format";
char *buff;
char testid[300];
i4 ise,ibe, next_order, next_tmp, bufsize, iorder;
pos hdrinit=2000, safety=100;
pos hdrsize, oldhdrsize, hdrsztmp, startloc;
int fotest;
/* Make certain (pos) is defined as 64bit int, otherwise, exit */

test_sizes();
fotest=0;


ibe = is_big_endian ();
ise = 0;
if (ibe == 0) ise = 1;    /* set flag for byteswapping if small endian */

/* Initially, open the file as read-only */

fp=fopen(fname,"rb");

if (fp == NULL) /* does file fname exist? */
{
	/* if not, print msg and exit */
        printf("sdf_insert: file %s does not exist; exiting\n",fname);
	fflush(stdout);
        exit(1);
}
else
{
	fotest=1;
	fgets(testid,strlen(sdfstr)+1,fp);
	/* Now, test to make sure the SDF file identifier string is there */
	/* If not, print msg and exit */
	if(strncmp(sdfstr,testid,strlen(sdfstr)))
	{
		fclose(fp);
		fotest=0;
		printf("sdf_insert: File %s is not an SDF file, exiting\n",
			fname);
		fflush(stdout);
		exit(1);
	}
}

if (fotest)
{
        /*close file if it is open */
        fclose(fp);
        fotest=0;
}

/* re-open file with both read/write permissions */

fp=fopen(fname,"rb+");

/* Now, read the locations for next writes of header info and data */

fseeko(fp, (pos) strlen(sdfstr)+1, SEEK_SET);
fread(&hdrpos,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrpos,(pos)1,sizeof(pos));

fread(&datapos,sizeof(pos),1,fp);
if(ise) byteswap((void *)&datapos,(pos)1,sizeof(pos));

fread(&next_order,sizeof(i4),1,fp);
if(ise) byteswap((void *)&next_order,(pos)1,sizeof(i4));

fread(&hdrsize,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrsize,(pos)1,sizeof(pos));

if(insert > next_order)
{
	printf("sdf_insert: can't insert before dataset %d in %s\n",
            insert,fname);
	fflush(stdout);
        fclose(fp);
        exit(1);
}

if( hdrpos >= (pos) (hdrsize - safety) )
{
/* now increase size of header area, and move the rest of the file up,
   adjust the current value of datapos accordingly */
 printf("sdf_insert: header area almost full, incrementing by hdrinit.\n");
 fflush(stdout);
 oldhdrsize=hdrsize;
 hdrsize+=hdrinit;
 hdrsztmp=hdrsize;
/*
 fseeko(fp,hdrpos-(pos)(sizeof(pos)),SEEK_SET);
*/
 fseeko(fp,(pos)(strlen(sdfstr)+1+2*sizeof(pos)+sizeof(i4)),SEEK_SET);
 if(ise) byteswap((void *)&hdrsztmp,(pos)1,sizeof(pos));
 fwrite(&hdrsztmp, sizeof(pos),1,fp);
 startloc=(pos) strlen(sdfstr)+1+sizeof(pos)+oldhdrsize;
 move_up(fp,hdrinit,startloc,datapos);
 fflush(fp);
/*
 fn=fileno(fp);
*/
 datapos += hdrinit;
/*
  ftruncate(fn, (pos) datapos);
*/
 sdf_file_truncate(fp, (pos) datapos);
}

bufsize=hdrsize-sizeof(pos)-sizeof(i4); 
buff=(char *)calloc(bufsize*sizeof(char),1);

off=(pos)0;
curpos = (pos) strlen(sdfstr)+1+sizeof(pos)+hdrsize;
for (iorder=0; iorder <= next_order; iorder++) 
/* loop must include next_order for case of insertion at eof */
{
	filepos=ftello(fp);
	if(iorder == insert) 
	{
		id_new->order=insert;
		size_ins = data_size(id_new);
                byte_ins = size_ins * id_new->nbpw;
		if(id_new->datatype == 'c')
		{
			byte_ins *= (pos)2;
		}
		/* Insert the new dataset, lengthen file, reposition */

		/* DEBUG
                printf("byte_ins = %d, curpos = %d, datapos = %d\n",
                    (i4)byte_ins,(i4)curpos,(i4)datapos);
		*/
		move_up(fp, byte_ins, curpos, datapos);
		fflush(fp);
/*
		fn=fileno(fp);
		ftruncate(fn, (pos) datapos + (pos) byte_ins );
*/
                sdf_file_truncate(fp, (pos) datapos + (pos) byte_ins );
		fseeko(fp,curpos,SEEK_SET);
                /* test for complex variable type */
		if (id_new->datatype != 'c') 
		{
			if(ise) byteswap((void *)data_new, 
                              (pos) size_ins, id_new->nbpw);
        		fwrite(data_new,id_new->nbpw,size_ins,fp); 
                        /* is this 64-bit compatible? */
			/* for non-destructive, should btye-swap back! */
			if(ise) byteswap((void *)data_new, 
                              (pos) size_ins, id_new->nbpw);
		}
		else
		{
			if(ise) byteswap((void *)data_new, 
                               (pos) 2*size_ins, id_new->nbpw);
		        fwrite(data_new,id_new->nbpw,(pos)2*size_ins,fp); 
                        /*kludge for complex var */
			/* for non-destructive, should btye-swap back! */
			if(ise) byteswap((void *)data_new, 
                               (pos) 2*size_ins, id_new->nbpw);
		}
		curpos=ftello(fp);  /* SHOULD THIS BE HERE? */
		/* insert new string for new dataset */
		fseeko(fp,filepos,SEEK_SET);
		data_to_buff(id_new,buff+off,&offset);
		off+=offset;

	}
	if(iorder < next_order)
	{
		id=(data_id *) string_to_data(fp);
	        datasize=data_size(id);
		if(iorder >= insert) /* right limit for this? */
		{
			id->order++;
		}

		data_to_buff(id,buff+off,&offset);
		off+=offset;
		if(id->datatype != 'c') /* test for complex var */
 		{
              		curpos += (pos) datasize*(id->nbpw);
        	}
        	else
        	{
			curpos += (pos) datasize*2*(id->nbpw);
			/* kludge for complex variables */
        	}

		free(id->label);
		free(id->dims);
		free(id);
	}

	/* make sure you free id_new in calling program */
}

hdrpos = strlen(sdfstr)+1 + 2*sizeof(pos) + sizeof(i4) + sizeof(pos);
fseeko(fp, hdrpos, SEEK_SET);
/* Gonna delete the \0 here */
/* fprintf(fp,"%s%c",buff,'\0'); */
fprintf(fp,"%s",buff);
hdrpos= (pos) ftello(fp);
datapos = curpos;
next_order++;

hdrtmp=hdrpos;
datatmp=datapos;
next_tmp=next_order;
hdrsztmp=hdrsize;

if(ise) byteswap((void *)&hdrtmp,(pos)1,sizeof(pos));
if(ise) byteswap((void *)&datatmp,(pos)1,sizeof(pos));
if(ise) byteswap((void *)&next_tmp,(pos)1,sizeof(i4));
if(ise) byteswap((void *)&hdrsztmp,(pos)1,sizeof(pos));

fseeko(fp, (pos) strlen(sdfstr)+1, SEEK_SET);

fwrite(&hdrtmp,sizeof(pos),1,fp);
fwrite(&datatmp,sizeof(pos),1,fp);
fwrite(&next_tmp,sizeof(i4),1,fp);
fwrite(&hdrsztmp,sizeof(pos),1,fp);

fclose(fp);
free(buff);

/* We're done */

return 0;
}

i4 sdf_replace(char *fname, i4 replace, data_id * id_new, void * data_new)
{

/* Replace dataset "replace" with descriptor id_new and data_new in file fname*/

FILE *fp;
data_id *id;
pos hdrpos, datapos, hdrtmp, datatmp, datasize, curpos, off, offset;
pos size_new, byte_new, byte_old, shiftup, shiftdn, filepos;
char *sdfstr="SDF format";
char *buff;
char testid[300];
i4 ise,ibe, next_order, next_tmp, bufsize, iorder;
pos hdrinit=2000, safety=100;
pos hdrsize, oldhdrsize, hdrsztmp, startloc;
int fotest;

/* Make certain (pos) is defined as 64bit int, otherwise, exit */

test_sizes();
fotest=0;

ibe = is_big_endian ();
ise = 0;
if (ibe == 0) ise = 1;    /* set flag for byteswapping if small endian */

/* Initially, open the file as read-only */

fp=fopen(fname,"rb");

if (fp == NULL) /* does file fname exist? */
{
	/* if not, print msg and exit */
        printf("sdf_replace: file %s does not exist; exiting\n",fname);
	fflush(stdout);
        exit(1);
}
else
{
	fotest=1;
	fgets(testid,strlen(sdfstr)+1,fp);
	/* Now, test to make sure the SDF file identifier string is there */
	/* If not, print msg and exit */
	if(strncmp(sdfstr,testid,strlen(sdfstr)))
	{
		fclose(fp);
		fotest=0;
		printf("sdf_replace: File %s is not an SDF file, exiting\n",
			fname);
		fflush(stdout);
		exit(1);
	}
}

if (fotest)
{
        /*close file if it is open */
        fclose(fp);
        fotest=0;
}

/* re-open file with both read/write permissions */

fp=fopen(fname,"rb+");

/* Now, read the locations for next writes of header info and data */

fseeko(fp, (pos) strlen(sdfstr)+1, SEEK_SET);
fread(&hdrpos,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrpos,(pos)1,sizeof(pos));

fread(&datapos,sizeof(pos),1,fp);
if(ise) byteswap((void *)&datapos,(pos)1,sizeof(pos));

fread(&next_order,sizeof(i4),1,fp);
if(ise) byteswap((void *)&next_order,(pos)1,sizeof(i4));

fread(&hdrsize,sizeof(pos),1,fp);
if(ise) byteswap((void *)&hdrsize,(pos)1,sizeof(pos));

if(replace >= next_order)
{
	printf("sdf_replace: can't replace dataset %d in %s\n",
            replace,fname);
	fflush(stdout);
        fclose(fp);
        exit(1);
}

if( hdrpos >= (pos) (hdrsize - safety) )
{
/* now increase size of header area, and move the rest of the file up,
   adjust the current value of datapos accordingly */
 printf("sdf_replace: header area almost full, incrementing by hdrinit.\n");
 fflush(stdout);
 oldhdrsize=hdrsize;
 hdrsize+=hdrinit;
 hdrsztmp=hdrsize;
/*
 fseeko(fp,hdrpos-(pos)(sizeof(pos)),SEEK_SET);
*/
 fseeko(fp,(pos)(strlen(sdfstr)+1+2*sizeof(pos)+sizeof(i4)),SEEK_SET);
 if(ise) byteswap((void *)&hdrsztmp,(pos)1,sizeof(pos));
 fwrite(&hdrsztmp, sizeof(pos),1,fp);
 startloc=(pos) strlen(sdfstr)+1+sizeof(pos)+oldhdrsize;
 move_up(fp,hdrinit,startloc,datapos);
 fflush(fp);
/*
 fn=fileno(fp);
*/
 datapos += hdrinit;
/*
  ftruncate(fn, (pos) datapos);
*/
 sdf_file_truncate(fp, (pos) datapos);
}

bufsize=hdrsize-sizeof(pos)-sizeof(i4); 
buff=(char *)calloc(bufsize*sizeof(char),1);

off=(pos)0;
curpos = (pos) strlen(sdfstr)+1+sizeof(pos)+hdrsize;
for (iorder=0; iorder < next_order; iorder++) 
{
	id=(data_id *) string_to_data(fp);
	filepos=ftello(fp);
        datasize=data_size(id); /* size of old dataset */
	if(iorder == replace) 
	{
		id_new->order=replace;
		size_new = data_size(id_new); /* size of new dataset */
                byte_new = size_new * id_new->nbpw;
		if(id_new->datatype == 'c')
		{
			byte_new *= (pos)2;
		}
		byte_old = datasize*id->nbpw;
		if(id->datatype == 'c')
		{
			byte_old*= (pos)2;
		}
		shiftup=byte_new - byte_old;
		if(shiftup < 0) 
		{
			/* if true, shrink file size */
			shiftdn = -shiftup;
			move_dn(fp,shiftdn,curpos,datapos);
			fflush(fp);
/*
			fn=fileno(fp);
			ftruncate(fn,(pos)datapos - (pos)shiftdn );
*/
			sdf_file_truncate(fp, (pos) datapos - (pos) shiftdn );
		}
		else if (shiftup > 0)
		{
			/* if true, expand file size */
			move_up(fp,shiftup,curpos,datapos);
			fflush(fp);
/*
			fn=fileno(fp);
			ftruncate(fn,(pos)datapos + (pos)shiftup );
*/
			sdf_file_truncate(fp, (pos) datapos + (pos) shiftup );
		}

		/* now must actually write the new data to the file */

		/* DEBUG
                printf("byte_ins = %d, curpos = %d, datapos = %d\n",
                    (i4)byte_ins,(i4)curpos,(i4)datapos);
		*/

		fseeko(fp,curpos,SEEK_SET); /* position file for write */

                /* test for complex variable type */
		if (id_new->datatype != 'c') 
		{
			if(ise) byteswap((void *)data_new, 
                              (pos) size_new, id_new->nbpw);
        		fwrite(data_new,id_new->nbpw,size_new,fp); 
                        /* is this 64-bit compatible? */
			/* for non-destructive, should btye-swap back! */
			if(ise) byteswap((void *)data_new, 
                              (pos) size_new, id_new->nbpw);
		}
		else
		{
			if(ise) byteswap((void *)data_new, 
                               (pos) 2*size_new, id_new->nbpw);
		        fwrite(data_new,id_new->nbpw,(pos)2*size_new,fp); 
                        /*kludge for complex var */
			/* for non-destructive, should btye-swap back! */
			if(ise) byteswap((void *)data_new, 
                               (pos) 2*size_new, id_new->nbpw);
		}
		curpos=ftello(fp);  /* SHOULD THIS BE HERE? */
		/* create new string for new dataset */
		fseeko(fp,filepos,SEEK_SET); 
		data_to_buff(id_new,buff+off,&offset);
		off+=offset;

	}

	else

	{

	        /* do this stuff only if dataset is not being replaced */
		data_to_buff(id,buff+off,&offset);
		off+=offset;
		if(id->datatype != 'c') /* test for complex var */
 		{
      			curpos += (pos) datasize*(id->nbpw);
       		}
       		else
       		{
			curpos += (pos) datasize*2*(id->nbpw);
			/* kludge for complex variables */
       		}
		/* end of stuff to be done if datset is not being replaced */

	}

	free(id->label);
	free(id->dims);
	free(id);

}
	/* remember to free id_new in calling program */

hdrpos = strlen(sdfstr)+1 + 2*sizeof(pos) + sizeof(i4) + sizeof(pos);
fseeko(fp, hdrpos, SEEK_SET);
/*  Gonna change this to get rid of \0 at end */
/* fprintf(fp,"%s%c",buff,'\0'); */
fprintf(fp,"%s",buff);
hdrpos= (pos) ftello(fp);
datapos = curpos;

hdrtmp=hdrpos;
datatmp=datapos;
next_tmp=next_order;
hdrsztmp=hdrsize;

if(ise) byteswap((void *)&hdrtmp,(pos)1,sizeof(pos));
if(ise) byteswap((void *)&datatmp,(pos)1,sizeof(pos));
if(ise) byteswap((void *)&next_tmp,(pos)1,sizeof(i4));
if(ise) byteswap((void *)&hdrsztmp,(pos)1,sizeof(pos));

fseeko(fp, (pos) strlen(sdfstr)+1, SEEK_SET);

fwrite(&hdrtmp,sizeof(pos),1,fp);
fwrite(&datatmp,sizeof(pos),1,fp);
fwrite(&next_tmp,sizeof(i4),1,fp);
fwrite(&hdrsztmp,sizeof(pos),1,fp);

fclose(fp);
free(buff);

/* We're done */

return 0;
}

void sdf_transpose_f77__(intf *indorderf77, intf *directionsf77, char * dtf77,
     intf *nbpwf77, intf *ndimf77, intf* dimsf77, void *dataf77, intf len_dtf77)
{
sdf_transpose_f77(indorderf77, directionsf77, dtf77, nbpwf77, ndimf77, dimsf77,
	dataf77, len_dtf77);
return;
}

void sdf_transpose_f77_(intf *indorderf77, intf *directionsf77, char * dtf77,
     intf *nbpwf77, intf *ndimf77, intf* dimsf77, void *dataf77, intf len_dtf77)
{
sdf_transpose_f77(indorderf77, directionsf77, dtf77, nbpwf77, ndimf77, dimsf77,
	dataf77, len_dtf77);
return;
}
void sdf_transpose_f77(intf *indorderf77, intf *directionsf77, char * dtf77,
     intf *nbpwf77, intf *ndimf77, intf* dimsf77, void *dataf77, intf len_dtf77)
{
char na[3]="NA";
data_id *id;
i4 *indorder;
i4 *directions;
i4 i,nlab;
id=malloc(sizeof(data_id));
id->ndim=(i4) *ndimf77;
id->nbpw=(i4) *nbpwf77;
id->datatype= (char) *dtf77;
nlab=strlen(na);
id->label=calloc(nlab+1,1);
strncpy(id->label,na,2);
id->dims=malloc(id->ndim*sizeof(pos));
indorder=malloc(id->ndim*sizeof(i4));
directions=malloc(id->ndim*sizeof(i4));

for (i=0;i<id->ndim;i++)
{
	*(id->dims+i)= (pos) *(dimsf77+id->ndim-1-i); /* flip dims for f77 */

	/*indorder starts at 1 in f77  - must flip order for f77. 
	 * Note also that definitions of indorder for C routine
         * are ndim-1- defs in fortran caller, +1 (because indorder starts at
         * 1 in fortran version).
         * I know this all seems baffling, but it gives correct results*/

	*(indorder+i)= (i4) id->ndim-1 - *(indorderf77+id->ndim-1-i)+1; 
	*(directions+i)= (i4) *(directionsf77+id->ndim-1-i); /*flip order f77 */
}

sdf_transpose(indorder, directions, id, dataf77);

/* must also change the dimsf77 array on output, since id->dims gets changed */
for (i=0; i < id->ndim; i++)
{
	*(dimsf77+i) = (intf) *(id->dims+id->ndim-1 -i); /* flip dims for f77 */
}

free(id->label);
free(id->dims);
free(id);
free(directions);
free(indorder);
return;
}

i4 sdf_transpose(i4 * indorder, i4 *directions, data_id *id, void *data)
/*

Use vacancy cycle tracking technique to convert multi-d array from original
index and dimension order to a permuted order, determined by the values 
pointed to by indorder.  The corresponding values of the directions
array determines whether or not the order of the new indices occurs
in their original order, or are reversed, depending on the sign of
the corresponding directions element.  After the function
is called, data will have been shuffled in place, and the order of id->dims
will be changed, to be consistent with indorder and the
sdf notation for array dimensions.
The vacancy cycle tracking technique of Ding (2001) 
[Chris H. Q. Ding, "An Optimal Index Reshuffle Algorithm for Multidimensional
Arrays and its Applications for Parallel Architectures", IEEE Transactions 
on Parallel and Distributed Systems, vol. 12, No. 3, pp 306-315, 2001]
allows the transformation
to be done in place, rather than copying from one array to a new array.

The speed of the reshuffle part of this procedure could be greatly speeded up
if we knew that the last index, or last set of indices, were the same as the
original.  In that case, the overall loop indexed by ir would be much smaller,
and the size of the memory transfers (the size of tmp) would be correspondingly
bigger.  The factor in each case would be the value of dim[last] (or for
several (q) dimensions together), the factor would be dim[last-q-1] x ... x
dim[last-2]*dim[last-1] - not sure that last bit is correct, but you get
the idea.  Roughing this out, which would involve testing indorder against
[0,...n-1] (ie newindorder, assuming it has no errors)
and then seeing how many contiguous last elements there were.
The dimension of the shuffle would also be reduced, so the values of ndim
passed into rindcalc and memcalcr would have to be changed as well.

Index reversal is done as a second step after
the re-shuffling has been accomplished.
It is possible the two steps could be combined into one, but that is too
complicated for me to contemplate.

*/
{
pos datasize,ir;
pos *indices_r, *indices;
/* pos memloc_c; */   /* not used here anymore */ 
/* pos memloc_r; */   /* not used here anymore */
pos istart,ilast,inext;
pos nbytes, nelem_cycle;
pos *newdims;
i4 *newindorder;
u4 *mask;
i4 i,reverse,badindex,trivial;
char *tmp;
test_sizes();
/* define the original order of dimensions and indices: */
datasize=data_size(id);
nbytes=(pos)id->nbpw;
if(id->datatype == 'c') nbytes*= (pos)2;
tmp=malloc(nbytes);
indices_r=malloc(id->ndim*sizeof(pos));
indices=malloc(id->ndim*sizeof(pos));
newindorder=malloc(id->ndim*sizeof(i4));
/* make sure the shuffled dimensions are correctly specified: */
for (i=0;i<id->ndim;i++)
{
	*(newindorder+i)=*(indorder+i);
			
}
/* sort copy of indorder array */
qsort((void *)newindorder,id->ndim,sizeof(i4),i4cmp); 
badindex=0;
trivial=1;
for (i=0;i<id->ndim;i++)
{
	if(*(newindorder+i) != (i4)i) 
	{
		badindex=1;
		printf("sdf_transpose: newindorder[%d] = %d\n",
			i,*(newindorder+i));
		fflush(stdout);
	}
	if(*(indorder+i) != (i4)i) trivial=0;
}

/* if the shuffled order of indorder and the original dim orders don't match,
   we're screwed, so bail and return without doing anything */

if(badindex)
{
	printf("sdf_transpose:  error in indorder, no transpose done\n");
	fflush(stdout);
	free(newindorder);
	free(tmp);
	free(indices_r);
	free(indices);
	return 1;
}

if(trivial == 0)

/* START OF RESHUFFLE: */

{
	newdims=malloc(id->ndim*sizeof(pos));
	mask=(u4 *)calloc(datasize/(sizeof(u4)*(pos)8)+(pos)1,sizeof(u4));
/* First, get array of new dimensions in the desired order: */
	for (i=0;i<id->ndim;i++)
	{
		newdims[i]=id->dims[indorder[i]];
	}
/* Loop over each gridpoint, doing the vacancy cycle if not done already */
	for (ir=0;ir<datasize;ir++)
	{
		while(readmask(mask,datasize,ir) == 0)
		{
			istart=ir;
			ilast=istart;
			rindcalc(id->ndim,newdims,indices_r,ir); 
			for (i=0;i<id->ndim;i++)
			{
				indices[indorder[i]]=indices_r[i];
			}
			memcalcr(id->dims,indices,id->ndim,&inext);
			memcpy((void *)tmp,(void *)((char *)data+inext*nbytes),
                               nbytes);
			nelem_cycle=(pos)0;
			writemask(mask,datasize,istart,1);
			while(inext != istart)
			{
				ilast=inext;
				rindcalc(id->ndim,newdims,indices_r,ilast); 
				for (i=0; i< id->ndim; i++)
				{
					indices[indorder[i]]=indices_r[i];
				}
				memcalcr(id->dims,indices,id->ndim,&inext);
				memcpy((void *)((char *)data+ilast*nbytes),
       		                         (void *)((char *)data+inext*nbytes), 
					   nbytes);
				nelem_cycle++;
				writemask(mask,datasize,ilast,1);
			}
                        /*
			 printf("nelem_cycle = %d\n",(i4)nelem_cycle);
			*/
			memcpy((void *)((char *)data+istart*nbytes),
				(void *)tmp,nbytes);
		}
	}
	
/* flip order of id->dims */
	for (i=0; i < id->ndim; i++)
	{
		*(id->dims+i)=(pos) *(newdims+i);
	}
	free(newdims);
	free(mask);
/*
END OF RESHUFFLE
*/
}

reverse=0;
for (i=0;i<id->ndim;i++)
{
	if(directions[i] < 0) reverse=1;
}
/*
Next block of code does the index reversal, if any of the elements of
the directions array are negative:

BEGIN INDEX REVERSAL
*/

if(reverse)
{
	/* re-allocate and reset bit mask: */
	mask=(u4 *)calloc(datasize/(sizeof(u4)*(pos)8)+(pos)1,sizeof(u4));
	for (ir=0;ir<datasize;ir++)
	{
		while(readmask(mask,datasize,ir) == 0)
		{
			rindcalc(id->ndim,id->dims,indices_r,ir);
			for (i=0;i < id->ndim; i++)
			{
				indices[i] = (directions[i] < 0) 
				? id->dims[i]-(pos)1-(pos)indices_r[i] 
				: indices_r[i];	
			}
			memcpy((void *)tmp,(void *)((char *)data+ir*nbytes),
				nbytes);
        		memcalcr(id->dims,indices,id->ndim,&inext);
                        memcpy((void *)((char *)data+ir*nbytes),
                                (void *)((char *)data+inext*nbytes), nbytes);
			memcpy((void *)((char *)data+inext*nbytes),
				(void *)tmp,nbytes);
			writemask(mask,datasize,ir,1);
			writemask(mask,datasize,inext,1);
		}
	}
	free(mask);
}
/* END INDEX REVERSAL */

/* clean up and return */
free(newindorder);
free(tmp);
free(indices);
free(indices_r);
return 0;
}

i4 rindcalc(i4 ndim, pos *dims, pos *indices, pos memloc_r)
{
/*  This function returns an array of indices corresponding to the memory
    location memloc_r (assuming row major) and dimensions dims.
    Assumes memory for pointers dims and indices are allocated in the
    caller.
*/
pos bigq;
i4 i;
bigq = memloc_r;
for (i=ndim-1;i>=0;i--)
{
	*(indices+i) = bigq % *(dims+i);
        bigq /= *(dims+i);
}
return 0;
}

i4 cindcalc(i4 ndim, pos *dims, pos *indices, pos memloc_c)
{
/*  This function returns an array of indices corresponding to the memory
    location memloc_c (assuming column major) and dimensions dims.
    Assumes memory for pointers dims and indices are allocated in the
    caller.
*/
pos bigm;
i4 i;
bigm = memloc_c;
for (i=0;i<ndim;i++)
{
	*(indices+i) = bigm % *(dims+i);
        bigm /= *(dims+i);
}
return 0;
}

i4 memcalcr(pos *dims, pos *indices, i4 ndim, pos *memloc_r)
{
/*
Given dimensions dims and indices indices, allocated in caller, this function
calculates the row major equivalent 1-d location in memory
(memloc_r) for a multidimensional array element.  Code to compute the
column major location memloc_c as well is here, but commented out to speed
things up.
*/
pos bigq; 
/*
pos bigm;
*/
i4 i;
bigq =(pos)0;
/*
 bigm =(pos)0;
*/
for (i=0;i<ndim;i++)
{
	bigq=*(dims+i)*bigq+*(indices+i);
	/*
         bigm=*(dims+ndim-1-i)*bigm+*(indices+ndim-1-i);
	*/
}
/*
 *memloc_c=bigm;
*/
*memloc_r=bigq;
return 0;
}

i4 memcalcc(pos *dims, pos *indices, i4 ndim, pos *memloc_c)
{
/*
Given dimensions dims and indices indices, allocated in caller, this function
calculates the column major equivalent 1-d location in memory
(memloc_c) for a multidimensional array element.  Row-major analogues are
included, but commented out.
*/
/*
pos bigq; 
*/
pos bigm;
i4 i;
/*
bigq =(pos)0;
*/
bigm =(pos)0;
for (i=0;i<ndim;i++)
{
	/*
	bigq=*(dims+i)*bigq+*(indices+i);
	*/
	bigm=*(dims+ndim-1-i)*bigm+*(indices+ndim-1-i);
}
*memloc_c=bigm;
/*
*memloc_r=bigq;
*/
return 0;
}

i4 readmask(u4 * mask, pos sizem, pos loc)

/* reads the single bit corresponding to location loc for a mask of
   size sizem.  Assumes user has allocated enough memory for mask in caller,
   and which is assumed to be divided into an array of 32 bit unsigned integers.
   Here is an example:

   mask=calloc(sizem/(sizeof(u4)*(pos)8)+(pos)1,sizeof(u4));

   */
{
i4 result;
pos intindex;
u4 bitindex, baseint;
if(loc >= sizem)
{
   printf("readmask: loc is out of range\n");
   exit(1);
}
intindex= (pos)(loc / (sizeof(u4)*(pos)8));
bitindex= (u4)(loc % (sizeof(u4)*(pos)8));
baseint=*(mask+intindex);
result= (baseint & ((u4)1 << bitindex)) ? 1 : 0;
return result;
}

i4 writemask(u4 *mask, pos sizem, pos loc, i4 value)
{
/* Write out a single bit according to value; if value=0, a 0 will be written,
   for any other value, a 1 will be written.  This is done for a bit at
   location loc, corresponding to a total bit mask length of sizem.  
   It is assumed
   the user has allocated enough space for mask in caller,
   which is assumed to be an array of unsigned 32 bit integers 
   Here is an example:

   mask=calloc(sizem/(sizeof(u4)*(pos)8)+(pos)1,sizeof(u4));

   */


pos intindex;
i4 result;
u4 baseint,bitindex;
if(loc >= sizem)
{
   printf("writemask: loc is out of range\n");
   exit(1);
}
intindex= (pos)(loc / (sizeof(u4)*(pos)8));
bitindex= (u4)(loc % (sizeof(u4)*(pos)8));
baseint=*(mask+intindex);
result = value ? (baseint | ((u4)1 << bitindex)) : 
        (baseint & ~((u4)1 << bitindex)); /* is there a simpler way? */
*(mask+intindex)=result;
return 0;

}
i4 byteswap (unsigned char *arr, pos arrsize, i4 nbpw)
/* Pretty simple:  arr is input array, which is byte-swapped in place,
   nbpw is the number of bytes per word, and arrsize is the size of the array
   (in units of nbpw bytes).  It is assumed that arr has
   already have been correctly defined and allocated in the calling program. */
{
  pos i, j;
  unsigned char temp;
/* declare temp a register variable? C newsgroups vary on wisdom of this */

  if(nbpw == 1) return 0; /*no point in byteswapping an array of bytes */

  for (i = (pos) 0; i < (pos) arrsize; i++) /* loop over the array elements */
    {
      for (j = (pos)0; j < (pos) nbpw/2; j++)/* loop over bytes in an element */
        {
          temp = *(arr + i*(pos) nbpw + ((pos) nbpw - j - (pos)1));
          *(arr + i*(pos) nbpw + ((pos)nbpw - j - (pos)1)) = 
                 *(arr +i*(pos)nbpw + j);
          *(arr + i* (pos)nbpw + j) = temp;
        }
    }
  return 0;
}

i4 is_big_endian ()
/* This function returns 1 if it is a big endian machine, 0 otherwise */
{
  unsigned char fakeword[4] = { 0xFF, 0x00, 0xFF, 0x00 };
  i4 realword = *((i4 *) fakeword);
  if (realword == 0xFF00FF00)
    {
      return 1;
    }
  else
    {
      return 0;
    }
}

void ***** sdf_1d_to_5d(data_id *id, void *data)
{

/* This function returns a 5d dynamically allocated array, based on the
   values of nbpw, ndim, and dims in the id structure.  The values of the
   array are set by the 1-d array data read in from sdf_read */

        char *****data5d;
	i4 sizeof_elem;
	i4 sizeof_ptr;
	pos i,j,k,l;
        sizeof_elem=id->nbpw; /* compute size of data element */
	if(id->datatype == 'c')
	{
		sizeof_elem*=2;
	}
	sizeof_ptr=sizeof(char *); /* compute size of pointer  */
	test_sizes();
	if(id->ndim != 5)
	{
		printf("sdf_mk_5d: id->ndim != 5; returning NULL\n");
		fflush(stdout);
		return NULL;
	}
        data5d=sdf_malloc(id->dims[0]*sizeof_ptr);
        for (i=0;i<id->dims[0];i++)
	{
		data5d[i]=malloc(id->dims[1]*sizeof_ptr);

		for (j=0;j<id->dims[1];j++)
		{
			data5d[i][j]=malloc(id->dims[2]*sizeof_ptr);
			for (k=0;k<id->dims[2];k++)
			{
				data5d[i][j][k]=malloc(id->dims[3]*sizeof_ptr);
			}
		}
	}

        data5d[0][0][0][0]=data;

        /* compute addresses for each start of each column */
        for (i=1;i<id->dims[0];i++)
	{
		data5d[i][0][0][0]=data5d[i-1][0][0][0]+id->dims[4]*id->dims[3]*
			id->dims[2]*id->dims[1]*sizeof_elem;
	}

	for (i=0;i<id->dims[0];i++)
	{
		for (j=1;j<id->dims[1];j++)
		{
			data5d[i][j][0][0]=data5d[i][j-1][0][0]+id->dims[4]*
				id->dims[3]*id->dims[2]*sizeof_elem;
		}
	}
	
	for (i=0;i<id->dims[0];i++)
	{
		for (j=0;j<id->dims[1];j++)
		{
			for (k=1;k<id->dims[2];k++)
			{
				data5d[i][j][k][0]=data5d[i][j][k-1][0] +
					id->dims[4]*id->dims[3]*sizeof_elem;
			}
		}
	}
	
	for (i=0;i<id->dims[0];i++)
	{
		for (j=0;j<id->dims[1];j++)
		{
			for (k=0;k<id->dims[3];k++)
			{
				for (l=1;l<id->dims[4];l++)
				{
					data5d[i][j][k][l]=data5d[i][j][k][l-1]
					+ id->dims[4]*sizeof_elem;
				}
			}
		}
	}

	return (void *****)data5d;
}
void **** sdf_1d_to_4d(data_id *id, void *data)
{

/* This function returns a 4d dynamically allocated array, using the
   array values from the 1d array data returned e.g. from sdf_read */

        char ****data4d;
	i4 sizeof_elem;
	i4 sizeof_ptr;
	pos i,j,k;
        sizeof_elem=id->nbpw; /* compute size of data element */
	if(id->datatype == 'c')
	{
		sizeof_elem*=2;
	}
	sizeof_ptr=sizeof(char *); /* compute size of pointer  */
	test_sizes();
	if(id->ndim != 4)
	{
		printf("sdf_1d_to_4d: id->ndim != 4; returning NULL\n");
		fflush(stdout);
		return NULL;
	}
        data4d=sdf_malloc(id->dims[0]*sizeof_ptr);
        for (i=0;i<id->dims[0];i++)
	{
		data4d[i]=malloc(id->dims[1]*sizeof_ptr);

		for (j=0;j<id->dims[1];j++)
		{
			data4d[i][j]=malloc(id->dims[2]*sizeof_ptr);
		}
	}

	data4d[0][0][0]=data;

        /* compute addresses for each start of each column */
        for (i=1;i<id->dims[0];i++)
	{
		data4d[i][0][0]=data4d[i-1][0][0]+id->dims[3]*id->dims[2]*
			id->dims[1]*sizeof_elem;
	}
	for (i=0;i<id->dims[0];i++)
	{
		for (j=1;j<id->dims[1];j++)
		{
			data4d[i][j][0]=data4d[i][j-1][0]+id->dims[3]*
				id->dims[2]*sizeof_elem;
		}
	}
	
	for (i=0;i<id->dims[0];i++)
	{
		for (j=0;j<id->dims[1];j++)
		{
			for (k=1;k<id->dims[2];k++)
			{
				data4d[i][j][k]=data4d[i][j][k-1] +
					id->dims[3]*sizeof_elem;
			}
		}
	}

	return (void ****)data4d;
}

void *** sdf_1d_to_3d(data_id *id, void *data)
{

/* This function returns a 3-d dynamically allocated array, whose
   elements also coincide in memory with the elements of the 1-d array 
   returned by a call to sdf_read */

        char ***data3d;
	i4 sizeof_elem;
	i4 sizeof_ptr;
	pos i,j;
        sizeof_elem=id->nbpw; /* compute size of data element */
	if(id->datatype == 'c')
	{
		sizeof_elem*=2;
	}
	sizeof_ptr=sizeof(char *); /* compute size of pointer  */
	test_sizes();
	if(id->ndim != 3)
	{
		printf("sdf_1d_to_3d: id->ndim != 3; returning NULL\n");
		fflush(stdout);
		return NULL;
	}
        data3d=sdf_malloc(id->dims[0]*sizeof_ptr);
        for (i=0;i<id->dims[0];i++)
	{
		data3d[i]=malloc(id->dims[1]*sizeof_ptr);
	}

        data3d[0][0]=data;  /* here is where the data from sdf_read is */

        /* compute addresses for each start of each column */
        for (i=1;i<id->dims[0];i++)
	{
		data3d[i][0]=data3d[i-1][0]+id->dims[2]*id->dims[1]*sizeof_elem;
	}
	for (i=0;i<id->dims[0];i++)
	{
		for (j=1;j<id->dims[1];j++)
		{
			data3d[i][j]=data3d[i][j-1]+id->dims[2]*sizeof_elem;
		}
	}
	return (void ***) data3d;
}

void ** sdf_1d_to_2d(data_id *id, void *data)
{

/* This function returns a 2-d dynamically allocated array, whose
   elements also coincide in memory with the elements of the 1-d array 
   returned by a call to sdf_read */

        char **data2d;
	i4 sizeof_elem;
	i4 sizeof_ptr;
	pos i;
        sizeof_elem=id->nbpw; /* compute size of data element */
	if(id->datatype == 'c')
	{
		sizeof_elem*=2;
	}
	sizeof_ptr=sizeof(char *); /* compute size of pointer  */
	test_sizes();
	if(id->ndim != 2)
	{
		printf("sdf_1d_to_2d: id->ndim != 2; returning NULL\n");
		fflush(stdout);
		return NULL;
	}
        data2d=sdf_malloc(id->dims[0]*sizeof_ptr);

        data2d[0]=data;  /* here is where the data from sdf_read is */

        /* compute addresses for each start of each column */
        for (i=1;i<id->dims[0];i++)
	{
		data2d[i]=data2d[i-1]+id->dims[1]*sizeof_elem;
	}
	return (void **) data2d;
}

void ** sdf_mk_2d(data_id *id)
{

/* This function returns a 2-d dynamically allocated array, whose
   elements have npbw, ndim, and dims from the id structure */

        char **data2d;
	i4 sizeof_elem;
	i4 sizeof_ptr;
	pos i;
        sizeof_elem=id->nbpw; /* compute size of data element */
	if(id->datatype == 'c')
	{
		sizeof_elem*=2;
	}
	sizeof_ptr=sizeof(char *); /* compute size of pointer  */
	test_sizes();
	if(id->ndim != 2)
	{
		printf("sdf_mk_2d: id->ndim != 2; returning NULL\n");
		fflush(stdout);
		return NULL;
	}
        data2d=sdf_malloc(id->dims[0]*sizeof_ptr);

	data2d[0]=sdf_malloc(id->dims[0]*id->dims[1]*sizeof_elem);

        /* compute addresses for each start of each column */
        for (i=1;i<id->dims[0];i++)
	{
		data2d[i]=data2d[i-1]+id->dims[1]*sizeof_elem;
	}
	return (void **) data2d;
}
void *** sdf_mk_3d(data_id *id)
{

/* This function returns a 3d dynamically allocated array, based on the
   values of nbpw, ndim, and dims in the id structure. */

        char ***data3d;
	i4 sizeof_elem;
	i4 sizeof_ptr;
	pos i,j;
        sizeof_elem=id->nbpw; /* compute size of data element */
	if(id->datatype == 'c')
	{
		sizeof_elem*=2;
	}
	sizeof_ptr=sizeof(char *); /* compute size of pointer  */
	test_sizes();
	if(id->ndim != 3)
	{
		printf("sdf_mk_3d: id->ndim != 3; returning NULL\n");
		fflush(stdout);
		return NULL;
	}
        data3d=sdf_malloc(id->dims[0]*sizeof_ptr);
        for (i=0;i<id->dims[0];i++)
	{
		data3d[i]=malloc(id->dims[1]*sizeof_ptr);
	}

        data3d[0][0]=sdf_malloc(id->dims[0]*id->dims[1]*id->dims[2]*
		sizeof_elem);

        /* compute addresses for each start of each column */
        for (i=1;i<id->dims[0];i++)
	{
		data3d[i][0]=data3d[i-1][0]+id->dims[2]*id->dims[1]*sizeof_elem;
	}
	for (i=0;i<id->dims[0];i++)
	{
		for (j=1;j<id->dims[1];j++)
		{
			data3d[i][j]=data3d[i][j-1]+id->dims[2]*sizeof_elem;
		}
	}
	return (void ***)data3d;
}

void **** sdf_mk_4d(data_id *id)
{

/* This function returns a 4d dynamically allocated array, based on the
   values of nbpw, ndim, and dims in the id structure. */

        char ****data4d;
	i4 sizeof_elem;
	i4 sizeof_ptr;
	pos i,j,k;
        sizeof_elem=id->nbpw; /* compute size of data element */
	if(id->datatype == 'c')
	{
		sizeof_elem*=2;
	}
	sizeof_ptr=sizeof(char *); /* compute size of pointer  */
	test_sizes();
	if(id->ndim != 4)
	{
		printf("sdf_mk_4d: id->ndim != 4; returning NULL\n");
		fflush(stdout);
		return NULL;
	}
        data4d=sdf_malloc(id->dims[0]*sizeof_ptr);
        for (i=0;i<id->dims[0];i++)
	{
		data4d[i]=malloc(id->dims[1]*sizeof_ptr);

		for (j=0;j<id->dims[1];j++)
		{
			data4d[i][j]=malloc(id->dims[2]*sizeof_ptr);
		}
	}

        data4d[0][0][0]=sdf_malloc(id->dims[0]*id->dims[1]*id->dims[2]*
		id->dims[3]*sizeof_elem);

        /* compute addresses for each start of each column */
        for (i=1;i<id->dims[0];i++)
	{
		data4d[i][0][0]=data4d[i-1][0][0]+id->dims[3]*id->dims[2]*
			id->dims[1]*sizeof_elem;
	}
	for (i=0;i<id->dims[0];i++)
	{
		for (j=1;j<id->dims[1];j++)
		{
			data4d[i][j][0]=data4d[i][j-1][0]+id->dims[3]*
				id->dims[2]*sizeof_elem;
		}
	}
	
	for (i=0;i<id->dims[0];i++)
	{
		for (j=0;j<id->dims[1];j++)
		{
			for (k=1;k<id->dims[2];k++)
			{
				data4d[i][j][k]=data4d[i][j][k-1] +
					id->dims[3]*sizeof_elem;
			}
		}
	}

	return (void ****)data4d;
}

void ***** sdf_mk_5d(data_id *id)
{

/* This function returns a 5d dynamically allocated array, based on the
   values of nbpw, ndim, and dims in the id structure. */

        char *****data5d;
	i4 sizeof_elem;
	i4 sizeof_ptr;
	pos i,j,k,l;
        sizeof_elem=id->nbpw; /* compute size of data element */
	if(id->datatype == 'c')
	{
		sizeof_elem*=2;
	}
	sizeof_ptr=sizeof(char *); /* compute size of pointer  */
	test_sizes();
	if(id->ndim != 5)
	{
		printf("sdf_mk_5d: id->ndim != 5; returning NULL\n");
		fflush(stdout);
		return NULL;
	}
        data5d=sdf_malloc(id->dims[0]*sizeof_ptr);
        for (i=0;i<id->dims[0];i++)
	{
		data5d[i]=malloc(id->dims[1]*sizeof_ptr);

		for (j=0;j<id->dims[1];j++)
		{
			data5d[i][j]=malloc(id->dims[2]*sizeof_ptr);
			for (k=0;k<id->dims[2];k++)
			{
				data5d[i][j][k]=malloc(id->dims[3]*sizeof_ptr);
			}
		}
	}

        data5d[0][0][0][0]=sdf_malloc(id->dims[0]*id->dims[1]*id->dims[2]*
		id->dims[3]*id->dims[4]*sizeof_elem);

        /* compute addresses for each start of each column */
        for (i=1;i<id->dims[0];i++)
	{
		data5d[i][0][0][0]=data5d[i-1][0][0][0]+id->dims[4]*id->dims[3]*
			id->dims[2]*id->dims[1]*sizeof_elem;
	}

	for (i=0;i<id->dims[0];i++)
	{
		for (j=1;j<id->dims[1];j++)
		{
			data5d[i][j][0][0]=data5d[i][j-1][0][0]+id->dims[4]*
				id->dims[3]*id->dims[2]*sizeof_elem;
		}
	}
	
	for (i=0;i<id->dims[0];i++)
	{
		for (j=0;j<id->dims[1];j++)
		{
			for (k=1;k<id->dims[2];k++)
			{
				data5d[i][j][k][0]=data5d[i][j][k-1][0] +
					id->dims[4]*id->dims[3]*sizeof_elem;
			}
		}
	}
	
	for (i=0;i<id->dims[0];i++)
	{
		for (j=0;j<id->dims[1];j++)
		{
			for (k=0;k<id->dims[2];k++)
			{
				for (l=1;l<id->dims[3];l++)
				{
					data5d[i][j][k][l]=data5d[i][j][k][l-1]
					+ id->dims[4]*sizeof_elem;
				}
			}
		}
	}

	return (void *****)data5d;
}

void sdf_free_5d(data_id *id, void ***** data5d)
{
	pos i,j,k;
	sdf_free(data5d[0][0][0][0]);
	for (i=0;i<id->dims[0];i++)
	{
		for (j=0;j<id->dims[1];j++)
		{
			for(k=0;k<id->dims[2];k++)
			{
				sdf_free(data5d[i][j][k]);
			}

			sdf_free(data5d[i][j]);
		}

		sdf_free(data5d[i]);

	}
	sdf_free(data5d);
}

void sdf_free_4d(data_id *id, void **** data4d)
{
	pos i,j;
	sdf_free(data4d[0][0][0]);
	for (i=0;i<id->dims[0];i++)
	{
		for (j=0;j<id->dims[1];j++)
		{
			sdf_free(data4d[i][j]);
		}

		sdf_free(data4d[i]);

	}
	sdf_free(data4d);
}

void    sdf_free_3d(data_id *id, void *** data3d)
{

/* This function frees a 3d dynamically allocated array.  It arguments are 
   the id structure used to create the array and the array itself. */

	pos i;
	sdf_free(data3d[0][0]);
        for (i=0;i<id->dims[0];i++) 
	{
		sdf_free(data3d[i]);
	}
        sdf_free(data3d);
	return;
}

void    sdf_free_2d(data_id *id, void ** data2d)
{

/* This function frees a 3d dynamically allocated array.  It arguments are 
   the id structure used to create the array and the array itself. */

	sdf_free(data2d[0]);
        sdf_free(data2d);
	return;
}

data_id * sdf_cp_id(data_id *id_0)
{

/* function that returns a copy of the id structure that is its argument */

	data_id *id;
	i4 i,nlab;
	id=sdf_malloc(sizeof(data_id));
	id->order=id_0->order;
	id->nbpw=id_0->nbpw;
	id->ndim=id_0->ndim;
	id->datatype=id_0->datatype;
	nlab=strlen(id_0->label)+1;
	id->label=malloc(nlab*sizeof(char));
	id->dims=malloc(id->ndim*sizeof(pos));
	for (i=0;i<nlab;i++)
	{
		*(id->label+i)=*(id_0->label+i);	
	}
	for (i=0;i<id->ndim;i++) 
	{
		id->dims[i]=id_0->dims[i];
	}
	return id;
}
void sdf_free_id(data_id *id)
{
/* free the members of id, then free id itself.  It uses sdf_free so it
 * can be called from outside of sdf. */

	sdf_free(id->label);
	sdf_free(id->dims);
	sdf_free(id);

	return;
}
data_id *sdf_create_id(i4 order, char *label, char dtype, i4 nbpw, 
		i4 ndim, pos *dims)

/* Returns pointer to a data_id type structure, whose members are assigned
 * to values from the calling argument list.  The caller program needs to
 * free the structure with sdf_free_id once they are done using it, to avoid
 * a memory leak. */

{
	data_id * idnew;
	i4 i, nchar;

	idnew=sdf_malloc(sizeof(data_id));
	idnew->order=order;
	nchar=strlen(label)+1;
	idnew->label=sdf_malloc((pos)nchar);
	for (i=0;i<nchar-1;i++)
	{
		*(idnew->label+i) = *(label+i);
	}
	*(idnew->label+nchar-1) = '\0';
	idnew->datatype=dtype;
	idnew->nbpw=nbpw;
	idnew->ndim=ndim;
	idnew->dims=sdf_malloc((pos)ndim*sizeof(pos));
	for (i=0;i<ndim;i++)
	{
		*(idnew->dims+i)=*(dims+i);
	}
	return idnew;
}

void *sdf_malloc(pos nmemb)
{
	/* purpose of sdf_malloc and sdf_free are to provide self-consistent
	 * versions of malloc and free which can be called from outside the
	 * sdf library (like in a calling program) to avoid inconsistencies
	 * between versions of malloc/free from within sdf and versions that
	 * might use different versions of these functions in libraries
	 * outside of sdf.
	 */

	void *ptr;
	ptr=malloc(nmemb);
	if(ptr == NULL)
	{
		printf("sdf_malloc: memory allocation failed\n");
		fflush(stdout);
		exit(1);
	}
	return ptr;
}
void sdf_free(void *ptr)
{
	free(ptr);
	return;
}

void fflush_f77__()
{
fflush_f77();
return;
}
void fflush_f77_()
{
fflush_f77();
return;
}

void fflush_f77()
{
	/* This function just provides a fortran callable version of
	 * fflush(stdout) */
	fflush(stdout);
	return;
}

int i4cmp(const void *vp, const void *vq)
{
/* comparison fn for 2 i4 integers in qsort. */ 

        const i4 *p = vp;
        const i4 *q = vq;
        int result;
	/*
         i4 diff = *p - *q;
	*/
        if(*p < *q) result=(i4)(-1);
        if(*p > *q) result=(i4)1;
        if(*p == *q) result=(i4)0;
        return result;
}

void test_sizes()
{

/*
test to make sure all defined variables types have correct size 
*/

int i2size;
int i4size;
int u4size;
int i8size;
int f4size;
int f8size;
int possize;
int intfsize;
i2size=sizeof(i2);
i4size=sizeof(i4);
u4size=sizeof(u4);
i8size=sizeof(i8);
f4size=sizeof(f4);
f8size=sizeof(f8);
possize=sizeof(pos);
intfsize=sizeof(intf);
if (i2size != 2)
{
	printf("test_sizes: sizeof(i2) = %d != 2\n",i2size);
	fflush(stdout);
        exit(1);
}
if (i4size != 4)
{
	printf("test_sizes: sizeof(i4) = %d != 4\n",i4size);
	fflush(stdout);
        exit(1);
}
if (u4size != 4)
{
	printf("test_sizes: sizeof(u4) = %d != 4\n",u4size);
	fflush(stdout);
        exit(1);
}
if (f4size != 4)
{
	printf("test_sizes: sizeof(f4) = %d != 4\n",f4size);
	fflush(stdout);
        exit(1);
}
if (i8size != 8)
{
	printf("test_sizes: sizeof(i8) = %d != 8\n",i8size);
	fflush(stdout);
        exit(1);
}
if (f8size != 8)
{
	printf("test_sizes: sizeof(f8) = %d != 8\n",f8size);
	fflush(stdout);
        exit(1);
}
if (possize != 8)
{
	printf("test_sizes: sizeof(pos) = %d != 8\n",possize);
	fflush(stdout);
        exit(1);
}
if (intfsize != sizeof(double)/2)
{
	printf("test_sizes: sizeof(intf) = %d != sizeof(double)/2\n",intfsize);
        printf("NOTE: above test is the fortran standard for the size of an\n");
        printf("integer.  If you use a different int size, you must change\n");
        printf("this test in test_sizes.\n");

        /* NOTE: above test is the fortran standard for the size of an
        integer.  If you use something else, you'll need to change this test
        accordingly! */

	fflush(stdout);
        exit(1);
}

}

#ifdef WIN32

/* define strtok_r in terms of strtok; 
   fseeko, ftello for MS windoze, using 
   fsetpos, fgetpos, which are
   supposedly large-file compatible, unlike fseek and ftell. */
char *strtok_r(char* s1, char* sep, char** last)
{
	/* supposedly strtok itself is threadsafe in Windoze, so just use it */
	return strtok(s1, sep);
}

int fseeko(FILE *fd, pos offset, int whence)
{
        /* note -- whence implicitly assumed = SEEK_SET
         * therefore, this isn't going to work to e.g. get filesize
         * where you want to set whence to SEEK_END 
	 * But this could be fixed if you invoke the _stati64 function
	 * to find filesize, then compute offset relative to SEEK_END
         */

        return fsetpos(fd, &offset);
}
pos ftello(FILE *fd)
{
        pos offset;
        fgetpos(fd, &offset);
        return offset;
}
int sdf_file_truncate(FILE *fd, pos size)
{
/* Windoze version of 64-bit compatible file_truncate */
	int fn;
        pos sizetest;
	HANDLE fh; /* a windoze file handle definition */
	fseeko(fd,size,SEEK_SET);
        sizetest=ftello(fd);
        if(sizetest != size)
	{
		char temp[300];
		printf("sdf_file_truncate win: ftello value != fseeko value\n");
 		output_int64(temp,sizetest);
		printf("sdf_file_truncate win: sizetest = %s\n",temp);
		fflush(stdout);
                fclose(fd);
		exit(1);
	}
	fn=fileno(fd);
	fflush(fd);
        fh= (HANDLE) _get_osfhandle(fn); /* get windoze file handle */
        if(SetEndOfFile(fh) ==0) /* this sets 64-bit compatible file size */
	{
		char buf[100];
		sprintf(buf, "SetEndOfFile failed in sdf_file_truncate: %d\n",
				GetLastError() );
                printf("%s\n",buf);
		fflush(stdout);
		return 1;
	}
	fflush(fd);
	return 0;
}
#else
int sdf_file_truncate(FILE *fd, pos size)
/* Unix-like version of sdf_file_truncate using ftruncate */
{
	int fn;
	fn=fileno(fd);
	fflush(fd);
	if(ftruncate(fn,size))
	{
		printf("problem ftruncating the file in sdf_file_truncate\n");
		fflush(stdout);
		fclose(fd);
		exit(1);
	}
	fflush(fd);
	return 0;
}
#endif