/* pg2d.c - Using PGPLOT to display simple FITS histograms from by hister2d   */
/*            And other assorted programs that I've written                   */
/*          Current caveats:                                                  */
/*            Flips along the y axis only, and does so in a bastardized way   */
/*            Bug in -f makes images unplottable ??                           */

#include <string.h>
#include <stdio.h>
#include <math.h>
#include <ctype.h>
#include "fitsio.h"
#include <cpgplot.h>

int minmax(float [], int, float []);
void setvp();
void colors(int, float, float);

int main(int argc, char *argv[])
{
	fitsfile *fptr;   /* FITS file pointer, defined in fitsio.h           */
	int status = 0;   /* CFITSIO status value MUST be initialized to zero!*/
	int i, bitpix, naxis, ctype, itf, high, flipflag=0, rangeflag=0;
	int labelflag=1, pixflag=0, xboxflag=0, yboxflag=0, numflag=1;
	int scaleflag=1;
	long naxes[2] = {1,1}, fpixel[2] = {1,1};
	float xrange[2], yrange[2], range[2], xbox[2], ybox[2], tr[6];
	float *pixels;
	char filename[FLEN_FILENAME], trange[FLEN_FILENAME];
	char timestring[24], ra[15], dec[15];

	if(argc <  4){
		fprintf(stderr,"USAGE: pg2d INPUT IMGXFER(0 to 2) COLOR(-6 to 6)\n"); 
		fprintf(stderr,"\tExtra options:\n\t-f flips y axis\n");
		fprintf(stderr,"\t-l supresses plot labels\n");
		fprintf(stderr,"\t-n supresses axis number labels\n");
		fprintf(stderr,"\t-p forces coordinates to be in pixel space\n");
		fprintf(stderr,"\t-s supresses color scale wedge\n");
		fprintf(stderr,"\t-r=NNNN,NNNN defines color scaling limits\n");
		fprintf(stderr,"\t-x=NNNN,NNNN defines x window, PIXEL COORDINATES\n");
		fprintf(stderr,"\t-y=NNNN,NNNN defines y window, PIXEL COORDINATES\n");
		return(-1);}
	if( (sscanf(argv[2],"%d",&itf)!=1) || itf < 0 || itf > 2){
		fprintf(stderr,"Invalid Image Transfer requested:\n");
		fprintf(stderr,"\t0 = Linear\n\t1 = Logarithmic\n\t2 = Square-root\n");
		fprintf(stderr,"\t\tYou're on crack so I'm picking linear.\n");
		itf=0;
	}
	if( (sscanf(argv[3],"%d",&ctype)!=1) || fabs(ctype) < 1 || fabs(ctype) > 6){
		fprintf(stderr,"Invalid Colormap Requested\n");
		fprintf(stderr,"\tNegative values invert the colormap\n");
		fprintf(stderr,"\t1 = Greyscale\n\t2 = Rainbow\n\t3 = Heat\n\t");
		fprintf(stderr,"4 = Cool\n\t5 = IRAF\n\t6 = AIPS\n");
		fprintf(stderr,"\t\tYou're on crack so I'm picking greyscale.\n");
		ctype=1;
	}
	strcpy(filename, argv[1]);

	/* Check to see if arguments read in correct formats                  */
	if(argc >= 4){
		for(i=4; i<argc; i++){
			// Flip flag
			if(argv[i][0]=='-' && 'f'==tolower(argv[i][1]) && 
			  argv[i][1]!='\0'){
				strcpy(filename,strncat(argv[1],"[*,-*]",6));
				flipflag=1;
			}
			// Label flag
			if(argv[i][0]=='-' && 'l'==tolower(argv[i][1]) && 
			  argv[i][1]!='\0'){
				labelflag=0;
			}
			// Pixel coords flag
			if(argv[i][0]=='-' && 'p'==tolower(argv[i][1]) && 
			  argv[i][1]!='\0'){
				pixflag=1;
			}
			// Supress number labels
			if(argv[i][0]=='-' && 'n'==tolower(argv[i][1]) && 
			  argv[i][1]!='\0'){
				numflag=0;
			}
			// Supress color scale wedge
			if(argv[i][0]=='-' && 's'==tolower(argv[i][1]) && 
			  argv[i][1]!='\0'){
				scaleflag=0;
			}
			// Range Constraints
			if(argv[i][0]=='-' && 'r'==tolower(argv[i][1]) && 
			  argv[i][1]!='\0'){
				strncpy(trange,argv[i]+3,strlen(argv[i])-2);
				if(sscanf(trange,"%f,%f",&range[0],&range[1])!=2){
					fprintf(stderr,"Nonfatal error: Malformed range\n");
					fprintf(stderr,"\t-r=NNNN,NNNN defines color scaling\n");
					rangeflag=0;
				}
				else{
					rangeflag=1;
				}
			}
			// Windowing function, in X pixel coordinates
			if(argv[i][0]=='-' && 'x'==tolower(argv[i][1]) && 
			  argv[i][1]!='\0'){
				strncpy(trange,argv[i]+3,strlen(argv[i])-2);
				if(sscanf(trange,"%f,%f",&xbox[0],&xbox[1])!=2){
					fprintf(stderr,"Nonfatal error: Malformed argument?\n");
					fprintf(stderr,"\t-x=NNNN,NNNN defines x axis range\n");
					xboxflag=0;
				}
				else{
					pixflag=1;
					xboxflag=1;
				}
			}
			// Windowing function, in X pixel coordinates
			if(argv[i][0]=='-' && 'y'==tolower(argv[i][1]) && 
			  argv[i][1]!='\0'){
				strncpy(trange,argv[i]+3,strlen(argv[i])-2);
				if(sscanf(trange,"%f,%f",&ybox[0],&ybox[1])!=2){
					fprintf(stderr,"Nonfatal error: Malformed argument?\n");
					fprintf(stderr,"\t-y=NNNN,NNNN defines x axis range\n");
					yboxflag=0;
				}
				else{
					pixflag=1;
					yboxflag=1;
				}
			}
		}
	}
	
	if (!fits_open_file(&fptr, filename, READONLY, &status)) {
		if (!fits_get_img_param(fptr, 2, &bitpix, &naxis, naxes, &status)) {
			if (naxis >= 3 || naxis == 0){
				printf("Fatal Error - Only 1D or 2D images are supported\n"); 
				return(-1);}
			else{
				pixels = (float *) malloc(naxes[0]*naxes[1] * sizeof(float));
				if (pixels == NULL) {
					printf("Fatal Error - Memory allocation error\n");
					return(-1);}
				fits_read_pix(fptr,TFLOAT,fpixel,naxes[0]*naxes[1],NULL,
								pixels,NULL,&status);

				if(rangeflag==0){
					high=minmax(pixels,naxes[0]*naxes[1],range);
					fprintf(stderr,"range[0]: %g\trange[1]: %g @ pixnum %d\n",
							range[0], range[1], high);
				}
	
				status=cpgopen("?");
//				printf("status: %d\n",status);
				status=0;
				setvp();
//				cpgpap(6.5,1.0);
				cpgpap(10.0,1.0);
				cpgsitf(itf);
				colors(ctype,1.,0.5);
				/* YOU MUST SET THE TR MATRIX AFTER THE COLOR FUNCTIONS		*/
				/* X_plotted = tr[0] + tr[1]*I + tr[2]*J 					*/
				/* Y_plotted = tr[3] + tr[4]*I + tr[5]*J 					*/
				if(fits_read_key(fptr,TFLOAT,"CRVAL1",&tr[0],NULL,&status)!=0){
					tr[0]=0.;}
				if(fits_read_key(fptr,TFLOAT,"CDELT1",&tr[1],NULL,&status)!=0){
					tr[1]=1.;}
				tr[2]=0.;  // FORCE NO ROTATION
				if(fits_read_key(fptr,TFLOAT,"CRVAL2",&tr[3],NULL,&status)!=0){
					tr[3]=0.;}
				tr[4]=0.;  // FORCE NO ROTATION
				if(fits_read_key(fptr,TFLOAT,"CDELT2",&tr[5],NULL,&status)!=0){
					tr[5]=1.;}
				status=0.;
				
				if(pixflag){
					tr[0]=0; tr[1]=1; tr[2]=0;
					tr[3]=0; tr[4]=0; tr[5]=1;
				}

				//	Ideally the [*,-*] syntax should adjust the WCS info for us
				//		...But we are far far far away from the ideal situation
				//		It didn't change the ymin, so I will now
				if(flipflag){
					tr[3]+=fabs(tr[5])*naxes[1];}
				
				xrange[0]=tr[0];
				xrange[1]=tr[0]+(naxes[0]*tr[1]);
				yrange[0]=tr[3];
				yrange[1]=tr[3]+(naxes[1]*tr[5]);
				// Correct?  I don't know!  I think so?
				tr[0]-=tr[1]/2.;		tr[3]-=tr[5]/2.;

				if(xboxflag){
					xrange[0]=xbox[0]; xrange[1]=xbox[1];
				}
				if(yboxflag){
					yrange[0]=ybox[0]; yrange[1]=ybox[1];
				}

//printf("xrange: %f %f\nyrange: %f %f\ntr: %f %f %f %f %f %f\n",xrange[0],xrange[1],yrange[0],yrange[1],tr[0],tr[1],tr[2],tr[3],tr[4],tr[5]);

				cpgenv(xrange[0],xrange[1],yrange[0],yrange[1],0,-1);
				cpgimag(pixels,naxes[0],naxes[1],1,naxes[0],1,naxes[1],
						range[0],range[1],tr);
				cpgsch(0.6);
				if(labelflag){
					if(numflag){
						cpgbox("bcnts",0.,0,"bcnts",0.,0);
					}
					else{
						cpgbox("bcts",0.,0,"bcts",0.,0);
					}
					if(scaleflag){
						cpgwedg("RI",2.0,4.5,range[0],range[1],"ADU");
					}
					cpgsch(1.0);
					if(fits_read_key(fptr,TSTRING,"UTCSTAMP",
						&timestring,NULL,&status)!=0){
							fprintf(stderr,"UTCSTAMP not found in header!\n");
							timestring[0]='\0';
					}
					if(!pixflag){
						if(fits_read_key(fptr,TSTRING,"RA",
							&ra,NULL,&status)!=0){
								fprintf(stderr,"RA not found in header!\n");
								ra[0]='\0';
						}
						if(fits_read_key(fptr,TSTRING,"DEC",
							&dec,NULL,&status)!=0){
								fprintf(stderr,"DEC not found in header!\n");
								dec[0]='\0';
            		    }
					}

					cpgsch(0.8);
					cpgmtxt("B", 2.8, 0.0, 0.0, "NMSU/MSFC LCROSS GBOC");
					cpgmtxt("B", 4.0, 0.0, 0.0, timestring);

					cpgmtxt("B", 2.8, 0.685, 0.0, "Apache Point Observatory");
					cpgmtxt("B", 4.0, 0.745, 0.0, "ARC 3.5m Telescope");

					cpgsch(1.0);
//					cpglab("","",timestring);
/*
					cpgmtxt("B", 2.8, 0.0, 0.0, "RA: ");
					cpgmtxt("B", 2.8, 0.07, 0.0, ra);
					cpgmtxt("B", 4.0, 0.0, 0.0, "DEC: ");
					cpgmtxt("B", 4.0, 0.07, 0.0, dec);

					cpgmtxt("B", 2.8, 0.65, 0.0, "Apache Point Observatory");
					cpgmtxt("B", 4.0, 0.65, 0.0, "ARC 3.5m Telescope");
//					cpgmtxt("B", 4.8, -0.06, 0.0, timestring);
//					cpgmtxt("B", 6.0, -0.06, 0.0, ra);

					cpgsch(1.0);
					cpglab("","",timestring);
*/
				}
				else{
					cpgbox("bcts",0.,0,"bcts",0.,0);
					cpgwedg("RI",2.0,4.5,range[0],range[1],"");
					cpgsch(1.0);
				}
//				cpgiden();
				cpgend();
		   	}
		}
		fits_close_file(fptr, &status);
		free(pixels);
	} 
	if (status) fits_report_error(stderr, status);
	return(status);
}

int minmax(float a[], int n, float val[]){
	int i=0, b=0;
	val[0]=0.;	val[1]=0.;

	for(i=0; i<n; i++){
		if(a[i]<=val[0]){val[0]=a[i];}
		if(a[i]>=val[1]){val[1]=a[i]; b=i;}
	}
	return(b+1);
}

// Following two functions stolen from PGPlot example files
// Specifically pgdemo4.f

void setvp(){
	float d=0., vpx1=0., vpx2=0., vpy1=0., vpy2=0.;

	cpgsvp(0.,1.,0.,1.);
	cpgqvp(1,&vpx1,&vpx2,&vpy1,&vpy2);
	
	d=vpx2-vpx1; 
	if(d>=(vpy2-vpy1)){
		d=(vpy2-vpy1)/40.;}
	else {d/=40.;}

	vpx1 = vpx1 + 5.0*d;	vpx2 = vpx2 - 2.0*d;
	vpy1 = vpy1 + 8.0*d;	vpy2 = vpy2 - 2.0*d;
	cpgvsiz(vpx1,vpx2,vpy1,vpy2);
}

void colors(int type, float contra, float bright){
	float GL[2]={0.0, 1.0};
	float GR[2]={0.0, 1.0};
	float GG[2]={0.0, 1.0};
	float GB[2]={0.0, 1.0};

	float RL[9]={-0.5, 0.0, 0.17, 0.33, 0.50, 0.67, 0.83, 1.0, 1.7};
	float RR[9]={ 0.0, 0.0,  0.0,  0.0,  0.6,  1.0,  1.0, 1.0, 1.0};
	float RG[9]={ 0.0, 0.0,  0.0,  1.0,  1.0,  1.0,  0.6, 0.0, 1.0};
	float RB[9]={ 0.0, 0.3,  0.8,  1.0,  0.3,  0.0,  0.0, 0.0, 1.0};

	float HL[5]={ 0.0, 0.2, 0.4, 0.6, 1.0};
	float HR[5]={ 0.0, 0.5, 1.0, 1.0, 1.0};
	float HG[5]={ 0.0, 0.0, 0.5, 1.0, 1.0};
	float HB[5]={ 0.0, 0.0, 0.0, 0.3, 1.0};

	float CL[5]={ 0.0, 0.2, 0.4, 0.6, 1.0};
	float CR[5]={ 0.0, 0.0, 0.0, 0.3, 1.0};
	float CG[5]={ 0.0, 0.0, 0.5, 1.0, 1.0};
	float CB[5]={ 0.0, 0.5, 1.0, 1.0, 1.0};

	float WL[10]={ 0.0, 0.5, 0.5, 0.7, 0.7, 0.85, 0.85, 0.95, 0.95, 1.0};
	float WR[10]={ 0.0, 1.0, 0.0, 0.0, 0.3,  0.8,  0.3,  1.0,  1.0, 1.0};
	float WG[10]={ 0.0, 0.5, 0.4, 1.0, 0.0,  0.0,  0.2,  0.7,  1.0, 1.0};
	float WB[10]={ 0.0, 0.0, 0.0, 0.0, 0.4,  1.0,  0.0,  0.0, 0.95, 1.0};

	float AL[20]={ 0.0, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3, 
				   0.4, 0.4, 0.5, 0.5, 0.6, 0.6, 0.7, 
				   0.7, 0.8, 0.8, 0.9, 0.9, 1.0};
	float AR[20]={ 0.0, 0.0, 0.3, 0.3, 0.5, 0.5, 0.0, 
				   0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
				   1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
	float AG[20]={ 0.0, 0.0, 0.3, 0.3, 0.0, 0.0, 0.0, 
				   0.0, 0.8, 0.8, 0.6, 0.6, 1.0, 1.0, 
				   1.0, 1.0, 0.8, 0.8, 0.0, 0.0};
	float AB[20]={ 0.0, 0.0, 0.3, 0.3, 0.7, 0.7, 0.7, 
				   0.7, 0.9, 0.9, 0.0, 0.0, 0.0, 0.0, 
				   0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

	// Negative of colortable == -contra, 1-bright
	if (type < 0) {
		contra*=-1.;
		bright=1.-bright;
		type*=-1;
	}
		
	// Greyscale, Rainbow, Heat, IRAF, AIPS
	if (type == 1) cpgctab(GL, GR, GG, GB, 2, contra, bright);
	else if (type == 2) cpgctab(RL, RR, RG, RB, 9, contra, bright);
	else if (type == 3) cpgctab(HL, HR, HG, HB, 5, contra, bright);
	else if (type == 4) cpgctab(CL, CR, CG, CB, 5, contra, bright);
	else if (type == 5)	cpgctab(WL, WR, WG, WB, 10, contra, bright);
	else if (type == 6) cpgctab(AL, AR, AG, AB, 20, contra, bright);
	else cpgctab(GL, GR, GG, GB, 2, contra, bright);
}