/*
 * National Center for SuperComputing Applications, University of Illinois
 *
 * This NCSA software product is public domain software.  Permission
 * is hereby granted to do whatever you like with it. Should you wish
 * to make a contribution towards the production of this software, please
 * send us your comments about your experience with the software,  why
 * you liked or disliked it, how you use it, and most importantly, how it
 * helps your work. We will receive your comments at softdev@ncsa.uiuc.edu.
 *
 * Please send bug reports to bugs@ncsa.uiuc.edu
 *
 * Author: Eng-Whatt Toh, National Center for Supercomputing Applications
 *         ewtoh@ncsa.uiuc.edu
 */
#if ! defined(lint) && ! defined(LINT)
static char rcs_id[] = "$Id: td_miriad.c,v 1.2 1993/07/15 19:49:12 gbourhis Exp $";
#endif

/*
 *	File:		td_miriad.c
 *	Contents:	Miriad file calls
 *
 */

/* #include "td_com.h" */
#include "gr_com.h"
#include <math.h>
#include	<sys/types.h>
#include	<sys/stat.h>
#include	<stdio.h>
#include	<malloc.h>
#include	<fcntl.h>


/* ------------------------------------------------------- */
/*
 *	Load specified MIRIAD file into memory jng  nov/13/90
 */
int
td_MiriadLoad(db) /* MY */
A_Data_t	*db;
{
	Boolean first=TRUE;
	int i,j,k,ret=0;
	long 	imagesize, rowsize, tempsiz;

	int nitems;
	unsigned char* hbuf;

	unsigned char* mirheadopen();
   int 				fmir; /* miriad image file handle */
	int				nbytes;
	char 				mirname[1024];
	float32 			dummy; /* = first float in image file (not used) */
	float				*rowdata;

 printf("\n:td_mirload:src [%s] [%s]\n", db->dataName, db->pathName);

/* initializeto dummy values  */

   db->dims[0] = db->dims[1] = db->dims[2] = 1; /* ie a dot */ 
   db->rank = 3; db->min = 1000.0; db->max = -1000.0;

#ifdef MIRIADFMT

	hbuf = (unsigned char*) mirheadopen(db->pathName,&nitems);
   if (hbuf==NULL) {
		 gr_TextMsgOut("cannot read Miriad header\n");
       return(-1);
     }
   mirgetitem (hbuf,nitems,"naxis", &(db->rank));

	if (hbuf == NULL) return(-1);

	else {
   mirgetitem (hbuf,nitems,"naxis1",&db->dims[0]);
   mirgetitem (hbuf,nitems,"naxis2",&db->dims[1]);
   mirgetitem (hbuf,nitems,"naxis3",&db->dims[2]);
   printf("td_mirload: rank=%d dims= %d %d %d\n",
							db->rank,db->dims[0],db->dims[1],db->dims[2]);

	imagesize  = db->dims[0] * db->dims[1] * db->dims[2] * sizeof(float32);
   tempsiz  = mirimagesize(db->pathName);

	if ( tempsiz != imagesize ) {
		beep();
      printf(" File error: image file siz=%d mismatch with header (=%d)!!\n",
						 tempsiz,imagesize);
		return (-1);
		}
     

   /* otherwise, ok, proceed to malloc image buffer, read in headers */

		db->data   = (float32 ***)
			td_Malloc3Dfloat32(db->dims[0],db->dims[1],db->dims[2]);
		
	   if (db->data==NULL) printf("cannot malloc db->data\n");	

      mirgetitem (hbuf,nitems,"ctype1",db->label[0]);
      mirgetitem (hbuf,nitems,"ctype2",db->label[1]);
      mirgetitem (hbuf,nitems,"ctype3",db->label[2]);

      mirgetitem (hbuf,nitems,"crval1",&(db->crval[0]));
      mirgetitem (hbuf,nitems,"crval2",&(db->crval[1]));
      mirgetitem (hbuf,nitems,"crval3",&(db->crval[2]));

      mirgetitem (hbuf,nitems,"crpix1",&(db->crpix[0]));
      mirgetitem (hbuf,nitems,"crpix2",&(db->crpix[1]));
      mirgetitem (hbuf,nitems,"crpix3",&(db->crpix[2]));

      mirgetitem (hbuf,nitems,"cdelt1",&(db->cdelt[0]));
      mirgetitem (hbuf,nitems,"cdelt2",&(db->cdelt[1]));
      mirgetitem (hbuf,nitems,"cdelt3",&(db->cdelt[2]));


/* **jng  OPEN MIRIAD image file **/

	sprintf(mirname,"%s/image",db->pathName);
	if (-1 == (fmir = open(mirname,O_RDONLY)))
	   {printf("Unable to open %s -exiting.\n",mirname);
	    return (NULL);
	   }

	read(fmir,&dummy,sizeof(float32));

	rowdata  = (float32 *)td_Malloc1D(1,db->dims[0],
			(unsigned int)sizeof(float32), "temp float32");
	if(rowdata==NULL) printf("cannot malloc rowdata!\n");

	rowsize = db->dims[0] * sizeof(float32); /* size of one row of data */

/* ******************************** */
		if (db->data != NULL) 
		{
			for (i=0;i<db->dims[2];i++) {
				for (j=0;j<db->dims[1];j++) {

					/* get a row of miriad data, then flip and store in db->data */
					if (rowsize  != (nbytes = read(fmir,rowdata,rowsize )))
	   				{ printf("\t%s\n Try to read %i bytes, but only %i read\n",
		   								mirname, rowsize, nbytes);
						    return (-1);
	 					  }

					for (k=0;k<db->dims[0];k++) {
						db->data[k][j][i] = rowdata[k];
						if (first){ first = FALSE;
										db->min = db->max = db->data[0][0][0];
							}
						else if (db->min > db->data[k][j][i])
										db->min = db->data[k][j][i];
						else if (db->max < db->data[k][j][i])
										db->max = db->data[k][j][i];
						}
					}
				}   
            db->range = db->max - db->min;
			printf("td_mir: max=%f min=%f\n",db->max,db->min);

			close(fmir); /*** close miriad file */
			td_Free(rowdata);

			if (ret == -1)
			{
				sprintf(msg,
				"***ERROR: Cannot read 3D image from %s.\n",
					db->pathName);
				gr_TextMsgOut(msg);
				gr_TextMsgOut
				("Are you sure this is a 3D 32-bit floating point dataset?\n");
				td_Free3d(db->data);
			}
			else
				db->format = MIRIAD;
		}
		else
		{
			td_Free3d(db->data);
			ret = -1;
			sprintf(msg,
			"***ERROR: Not enough memory to load MIRIAD image %s.\n",
					db->pathName);
			gr_TextMsgOut(msg);
			gr_TextMsgOut
			("I suggest you quit NCSA X DataSlice and retry the loading.\n");
		}

		return(ret);
	}

#else
	return(-1);
#endif

}

/* ------------------------------------------------------- */

/* return Miriad header attributes as text */
/* calling routine must deallocate string when done */

char
*td_MiriadgetStats(db,nColors,nSplitColors)
A_Data_t    *db;
int         nColors,nSplitColors;
{
    int  i,j,ret;
    char *strng=NULL;

	float rtod = 180.0 / PI;
	float rtoh = 12.0  / PI;

   char s1[20],s2[20],s3[20]; /* temp strings */
	char tempstr[256];

#ifdef MIRIADFMT

    if (db->format != MIRIAD) return (NULL);
	 printf("in TDMIRIADSTATS\n");
    strng = td_Malloc1D(2500,1,sizeof(char),"td_Malloc1D:string");
    sprintf(strng,"Miriad File: %s\n", db->pathName);

		for (i=0;i<3;i++) {
            sprintf(tempstr,"%s: dim=%ld,\tlabel=%s\n",
                axesLabels[i],db->dims[i],db->label[i]);
			   strcat(strng, tempstr);
        }
    	sprintf(tempstr,"Max= %8.5f, Min=%8.5f\n", db->max,db->min);
		strcat(strng, tempstr);
    	sprintf(tempstr,
 			  "Header Data:\nNo.  crval        crpix        cdelt\n", strng);
		strcat(strng, tempstr);

		/* ---  Print Miriad header data, with formatting  -- */

 			/* --- first set of values (ra) */
           if( db->label[0][0]=='R' && db->label[0][1]=='A') 
 		        degtoang ((float) (db->crval[0]*rtoh),s1);
 		     else  
              degtoang ((float) (db->crval[0]*rtod),s1);
			  sprintf(s2,"%8.3f",(float) (db->crpix[0]));
 		     degtoang ((float) (db->cdelt[0]*rtod),s3);
           sprintf(tempstr," 1  %s  %s  %s\n",  s1, s2, s3);
			  strcat(strng, tempstr);

 			/* --- second set of values (dec) */
 			  degtoang ((float) (db->crval[1]*rtod),s1);
			  sprintf(s2,"%8.3f",(float) (db->crpix[1]));
 			  degtoang ((float) (db->cdelt[1]*rtod),s3);
     	     sprintf(tempstr," 2  %s  %s  %s\n",  s1, s2, s3);
			  strcat(strng, tempstr);

 			/* --- third set of values (vel) */
     	     sprintf(tempstr," 3  %8.3f  %8.3f  %8.3f\n", 
							db->crval[2], db->crpix[2], db->cdelt[2]);
			  strcat(strng, tempstr);

    db->range = db->max - db->min;
    db->rangeFrac = (float32)((float32)nColors/db->range);
    db->rangeFracSplit = (float32)((float32)nSplitColors/db->range);

#endif

    return(strng);
}

/* ------------------------------------------------------------------ */

/* returns a ptr to mem with header, and nitems */
unsigned char* mirheadopen(mirname,nitems)  /* newmiriad */
char *mirname;
int * nitems;

{
	struct	stat	mystat;
	char		fname[1024];
	unsigned char *hbuf;
	int		n, nbytes, fd;
	long val;
	char itemstr[90];
   char msg[280];

	strcpy(fname,mirname);
	strcat(fname,"/header");

	if (0 != stat(fname,&mystat))
	   {printf("Stat failed on %s.\n",fname);
		 gr_TextMsgOut("Not a miriad file!!\n"); beep();
		
	    return (NULL);
	   }

	if (0 == (n = mystat.st_size))
	   {printf("%s has zero-sized header - exiting.\n",fname);
	    return (NULL);
	   }

	if (0 == (hbuf = (unsigned char*) malloc(n)))
	   {printf("Unable to malloc %i bytes -exiting.\n",n);
	    return (NULL);
	   }
	
	if (-1 == (fd = open(fname,O_RDONLY)))
	   {printf("Unable to open %s -exiting.\n",fname);
	    return (NULL);
	   }

	if (n != (nbytes = read(fd,hbuf,n)))
	   {printf("File\t%s\nStat says %i bytes, %i read - exiting.\n",
		   fname, n, nbytes);
	    return (NULL);
	   }

	if (0 != close(fd))
	   {printf("Problem closing %s - continuing anyway.\n",fname);
	   }
	*nitems = 1+nbytes/32;

	return (hbuf);
}

/* ---------------------------------------------------------- */
/* returns value of itemstr in val. returns -1 if error. */
int mirgetitem (hbuf,nitems,itemstr,val) /* newmiriad */
unsigned char  hbuf[][32], *val;
int nitems;
char * itemstr;

 {
   int i,j, code, m;

	for(i=0;i<nitems ;i++) 
	  if (!strcmp(itemstr,(char *)hbuf[i]))  {
			code = (int) hbuf[i][19];
 printf("item=[%s] code=%d\n", itemstr, code);
			switch(code) {		
			case 1: /* string */
 				j=0;
				m = (int) ( hbuf[i][15] - 6);
				while (m-- >= -1) { val[j] = hbuf[i][20+j]; j++; }
				val[j] = '\0';
			break;

			case 2: /* int4 */	
				for(j=0;j<4;j++)  val[j] = (unsigned char) hbuf[i][20+j];
			break;

			case 3: /* int2 */	
				for(j=0;j<2;j++)  val[j] = (unsigned char) hbuf[i][20+j];
			break;

			case 4: /* float */	
				for(j=0;j<sizeof(float);j++)  val[j] = hbuf[i][20+j];
			break;

			case 5: /* double */	
				for(j=0;j<sizeof(double);j++)  val[j] = hbuf[i][24+j];
			break;

			default: printf("rubbish\n");	
   					return(-1); 
			break;
			} /* switch */
			return(0);
		}
	printf("no such item [%s]\n",itemstr);
   return(-1); 
}
/* ------------------------------------------------------- */
/* returns the miriad image file's byte size as found using stat() */
int mirimagesize(mirname) /* newmiriad */
char* mirname; {
	struct	stat	mystat;
	char		fname[1024];
	int		n, nbytes, fd;

	strcpy(fname,mirname);
	strcat(fname,"/image");

	if (0 != stat(fname,&mystat))
	   {printf("Stat failed on %s - exiting.\n",fname);
	    return (-1);
	   }

	if (0 == (n = mystat.st_size))
	   {printf("%s has zero-sized image - exiting.\n",fname);
	    return (-1);
	   }

	return(n-4); /* first 4 bytes in "/image" are not data at all */

}/* mirimagesize */
/* ------------------------------------------------------------------ */
/* converts degrees into string */
/* input: angle in degrees (float) */

degtoang (ang,ss) float ang; char ss[];  {
  char sign;
  int deg, min;
  float sec;
  float v;

  if(ang > 0.0) sign = '+'; else { sign='-';  ang = - ang; }

  v = ang;

  deg = (int) v;
  min = (v-(float)deg)*60;
  
  sec = (v - (float)deg - (float)min/60.0) *3600.0;
  sprintf(ss,"%c%d:%d:%05.2f",sign,deg,min,sec);
}

/* ------------------------------------------------------------------ */

/* jng 13-MAR-91. Used by contour hardcopy routine only */
/* if Miriad file, returns the x and y axes titles */
#ifdef MIRIADFMT
void
gr_getMiriadAxesTitles (dsplWin, xtitle, ytitle) 
A_DsplWind_t	*dsplWin;
char				* xtitle, *ytitle;

{
	A_Axes_t	 orient;
	A_Data_t	 *db=dsplWin->parent->parent->data;

	if (db->format != MIRIAD) return;

	orient=dsplWin->axesOrient;

	if (     ( orient.row  == 2 || orient.row  == -2 ) 
     		&& ( orient.col  == 1 || orient.col  == -1 ) 
     		&& ( orient.axis == 3 || orient.axis == -3 )  ) { /* XY  plane */
			strcpy(xtitle, db->label[0]);
			strcpy(ytitle, db->label[1]);
		}
	else 
	if (     ( orient.row  == 3 || orient.row  == -3 ) 
     		&& ( orient.col  == 1 || orient.col  == -1 ) 
     		&& ( orient.axis == 2 || orient.axis == -2 )  ) { /* XZ  plane */
			strcpy(xtitle, db->label[0]);
			strcpy(ytitle, db->label[2]);
		}
	else 
	if (     ( orient.row  == 3 || orient.row  == -3 ) 
     		&& ( orient.col  == 2 || orient.col  == -2 ) 
     		&& ( orient.axis == 1 || orient.axis == -1 )  ) { /* ZY  plane */
			strcpy(xtitle, db->label[1]);
			strcpy(ytitle, db->label[2]);
		}
}
#endif

/* ------------------------------------------------------------------ */
/* jng 13-MAR-91. Used by contour hardcopy routine only */
/* if Miriad file, returns the tick labels along the x and y axes */
/* FOR NOW, only 1 value (center) per axis returned */

#ifdef MIRIADFMT
void
gr_getMiriadAxesTickValues (dsplWin, xtick, ytick)
A_DsplWind_t	*dsplWin;
char 				*xtick, *ytick;

{
	A_Axes_t	 orient;
	A_Data_t	 *db=dsplWin->parent->parent->data;
	char		dummy[50];
	int selX, selY;

	if (db->format != MIRIAD) return;

	orient=dsplWin->axesOrient;

	if (     ( orient.row  == 2 || orient.row  == -2 ) 
     		&& ( orient.col  == 1 || orient.col  == -1 ) 
     		&& ( orient.axis == 3 || orient.axis == -3 )  ) { /* XY  plane */

			selX =  db->dims[0]/2-1; selY =  db->dims[1]/2-1;
			computeMiriadValues (selX, selY, db, orient, xtick, ytick, dummy);
		}
	else 
	if (     ( orient.row  == 3 || orient.row  == -3 ) 
     		&& ( orient.col  == 1 || orient.col  == -1 ) 
     		&& ( orient.axis == 2 || orient.axis == -2 )  ) { /* XZ  plane */

			selX =  db->dims[0]/2-1; selY =  db->dims[2]/2-1;
			computeMiriadValues (selX, selY, db, orient, xtick, ytick, dummy);
		}
	else 
	if (     ( orient.row  == 3 || orient.row  == -3 ) 
     		&& ( orient.col  == 2 || orient.col  == -2 ) 
     		&& ( orient.axis == 1 || orient.axis == -1 )  ) { /* ZY  plane */

			selX =  db->dims[1]/2-1; selY =  db->dims[2]/2-1;
			computeMiriadValues (selX, selY, db, orient, xtick, ytick, dummy);
		}
}
#endif
/* ------------------------------------------------------------------ */

/* jng 24-MAY-91. Used by contour hardcopy routine only */
/* if Miriad file, returns the axes hi and low end tick labels */

#ifdef MIRIADFMT
void
gr_getMiriadEndAxesTickValues (dsplWin, xtick0, ytick0, xtick1, ytick1)
A_DsplWind_t	*dsplWin;
char 				*xtick0, *ytick0;  /* low-end axes  value */
char 				*xtick1, *ytick1;  /* high-end axes value */

{
	A_Axes_t	 orient;
	A_Data_t	 *db=dsplWin->parent->parent->data;
	char		dummy[50];
	int selX, selY;

	if (db->format != MIRIAD) return;

	orient=dsplWin->axesOrient;

	selX =  selY = 0;
	computeMiriadValues (selX, selY, db, orient, xtick0, ytick0, dummy);

	if (     ( orient.row  == 2 || orient.row  == -2 ) 
     		&& ( orient.col  == 1 || orient.col  == -1 ) 
     		&& ( orient.axis == 3 || orient.axis == -3 )  ) { /* XY  plane */

			selX =  db->dims[0]; selY =  db->dims[1];
			computeMiriadValues (selX, selY, db, orient, xtick1, ytick1, dummy);
		}
	else 
	if (     ( orient.row  == 3 || orient.row  == -3 ) 
     		&& ( orient.col  == 1 || orient.col  == -1 ) 
     		&& ( orient.axis == 2 || orient.axis == -2 )  ) { /* XZ  plane */

			selX =  db->dims[0]; selY =  db->dims[2];
			computeMiriadValues (selX, selY, db, orient, xtick1, ytick1, dummy);
		}
	else 
	if (     ( orient.row  == 3 || orient.row  == -3 ) 
     		&& ( orient.col  == 2 || orient.col  == -2 ) 
     		&& ( orient.axis == 1 || orient.axis == -1 )  ) { /* ZY  plane */

			selX =  db->dims[1]; selY =  db->dims[2];
			computeMiriadValues (selX, selY, db, orient, xtick1, ytick1, dummy);
		}
}
#endif

/* ================================================================== */

#ifdef MIRIADFMT
/* display MIRIAD RA, DEC and VEL numerical values on spreadsheet window */
int
gr_DsplMiriadValues (dsplWin) 			/* jng 13-MAR-91 */
A_DsplWind_t	*dsplWin;
{
	A_Data_t	 *db=dsplWin->parent->parent->data;
	char 	 RAmsg[80], DECmsg[80], VELmsg[80];

	if (db->format != MIRIAD) return;

	/* display RA value in first box */
	/* display DEC value in second box */
	/* display VEL value in third box */

	printf("SELX: %d SELY: %d\n", dsplWin->selX, dsplWin->selY);

		computeMiriadValues (dsplWin->selX, dsplWin->selY, db, 
						&(dsplWin->axesOrient), RAmsg, DECmsg, VELmsg);

		gr_DialogSetValue(dsplWin->xDialog,RAmsg); 
		gr_DialogSetValue(dsplWin->yDialog,DECmsg); 
		gr_DialogSetValue(dsplWin->zDialog,VELmsg); 

}
#endif

/* ------------------------------------------------------------------ */
/*
	jng 8-apr-91
  Given a point (selX,selY)  on the display window,
  compute its corresponding Miriad values, return them as 3 strings.
  From the orient structure, it knows what the selX and selY corresponds to.

  Corrected to handle crdelts, crvals and crpix values properly.
 */

#ifdef MIRIADFMT
computeMiriadValues (selX, selY, db, orient, RAmsg, DECmsg, VELmsg) 

int			selX, selY; 	/* selected pt on display window */
A_Axes_t		*orient; 		/* orientation of that display window */
A_Data_t	 	*db;				/* miriad dataset */
char 			*RAmsg, *VELmsg, *DECmsg;  /* strings returned here */

{
	double		alpha,delta,vel,i,j,k;
	double		aa,dd,vv;
	double	tempcdelt; /* cosine adjustment for RA jng 5 june 91*/
	float rtod = 180.0 / PI;
	float rtoh = 12.0  / PI;

	if (db->format != MIRIAD) return;

	if (    ( orient->row  == 2 || orient->row  == -2 ) 
     		&& ( orient->col  == 1 || orient->col  == -1 ) 
     		&& ( orient->axis == 3 || orient->axis == -3 )  ) { /* XY  plane */
				aa = (double) selX;
				dd = (double) selY;
				vv = (double) orient->plane;
		}
	else 
	if (    ( orient->row  == 3 || orient->row  == -3 ) 
     		&& ( orient->col  == 1 || orient->col  == -1 ) 
     		&& ( orient->axis == 2 || orient->axis == -2 )  ) { /* XZ  plane */
				aa = (double) selX;
				vv = (double) selY;
				dd = (double) orient->plane;
		}
	else 
	if (    ( orient->row  == 3 || orient->row  == -3 ) 
     		&& ( orient->col  == 2 || orient->col  == -2 ) 
     		&& ( orient->axis == 1 || orient->axis == -1 )  ) { /* ZY  plane */
				dd = (double) selX;
				vv = (double) selY;
				aa = (double) orient->plane;
		}

   /*  - - - - - compute miriad header values - - - - - */

	vel   = db->crval[2] + db->cdelt[2] * (vv - db->crpix[2]);
	delta = db->crval[1] + db->cdelt[1] * (dd - db->crpix[1]);

	   if (db->label[0][0] =='R' && db->label[0][1]=='A') { 
			tempcdelt = db->cdelt[0]/cos(db->crval[1]);
			alpha = db->crval[0] + tempcdelt * (aa - db->crpix[0]);
			}
		else
			alpha = db->crval[0] + db->cdelt[0] * (aa - db->crpix[0]);

  printf ("computeM: alpha=%f delta=%f vel=%f\n",alpha,delta,vel);

	/* convert MIRIAD values into RAmsg, DECmsg and VELmsg */
	   if (db->label[0][0] =='R' && db->label[0][1]=='A') 
     		degtoang (alpha * rtoh,  RAmsg);
  		else
  		   degtoang (alpha * rtod, RAmsg);

	   degtoang (delta* rtod, DECmsg);
	   sprintf(VELmsg,"%f12.2",vel);
}
#endif

/* ================================================================== */

