/**************************************************************************
 *                                                                        *
 * Copyright (c) 1996 Michael Richmond and Richard Treffers               *
 *                                                                        *
 *                    This software may be copied and distributed for     *
 *                    educational, research and not for profit services   *
 *                    provided that this copyright and statement are      *
 *                    included in all such copies.                        *
 *                                                                        *
 **************************************************************************/


/*	generates phony target data
		PSTAR 
	
	12/11/92 - added field capability and changed nophoton to photon
			 - also stars are now specified by their peak value
			   instead of the integrated counts.
	12/16/92 - add evaluate to some inputs
	8/4/93  - add grid option
	8/24/93 - fix bug in photon gain
	3/3/94  - prototypes adjusted
	4/20/94 - change RAND_MAX
	1/8/96  - add 'help' option.  MWR
	1/29/96 - set a limit beyond which we don't evaluate contribution
	                    a star to the general background.  Will speed
	                    up execution a LOT.  MWR
*/

#include "pcvista.h"
#include "fits.h"
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <unistd.h>

#include <stdlib.h>	/* to get RAND_MAX */

#ifndef RAND_MAX
#include <limits.h>			/* needed to get INT_MAX */
#define RAND_MAX (INT_MAX)  /* watch out this may be wrong */
#endif

#ifndef M_PI
#define M_PI 3.14159265359
#endif

#define NBUF 80
#define MAXSTAR 1000 
#define DEFAULT_SIZE 64
#define DEFAULT_NAME "pstar"
#define NWIDTHS      4.0      /* don't bother to calculate the contribution */
                              /*   of a star to some pixel if that pixel is */
                              /*   this many FWHM from the star's center */

typedef int (*PFI)();

struct star_type{
	double height;
	double row;
	double col;
	double fwhm;
};

struct field_type{
	int nstar;
	double max; 
	double min;
	double fwhm;
};


#ifdef PROTO
static double gasdev(int *);
static double ran1(int *);
int main(int,char **);
static double gammln(double);
static int make_grid(struct star_type *s, struct field_type *f, 
	int nrow,int ncol);
static int make_field(struct star_type *s, struct field_type *f, 
	int nrow,int ncol);
static double poidev(double, int *);
static void show_options(FILE *fp);
int sort_by_row(struct star_type *a,struct star_type *b);
static double dgauss(double,double,double);
#else
static double gasdev();
static double ran1();
int main();
static double poidev();
static double gammln();
static int make_grid();
static int make_field();
static int sort_by_row();
static void show_options();
static double dgauss();
#endif

int
main(argc,argv)
int argc;
char *argv[];
{
	int i,j,k,warned=0;
	int photon=0;		/* whether to generate photon noise */
	int verbose_flag = 0;
	short data[NMAX];
	int nrow=DEFAULT_SIZE, ncol=DEFAULT_SIZE,idum,nstar=0,nfield;
	char *gotit,fname[NBUF],buf[NBUF],buf2[NBUF];
	FITS_HANDLE fh;
	double rms=0.,sky=0.,val;
	double rowgrad=0.,colgrad=0.,start;
	double gain=1.0;
	double dx, dy;
	double maxdist = 0.0, negmaxdist = 0.0, maxdist2 = 0.0;
	int got_field=0;
	int got_grid=0;
	struct star_type s[MAXSTAR];
	struct field_type f;

	if (argc == 1) {
		show_options(stderr);
		exit(1);
	}
	if (strcmp(argv[1], "help") == 0) {
		show_options(stdout);
		exit(0);
	}

	strcpy(fname,DEFAULT_NAME);

	for(i=1; i<argc; i++){
		if (keyword("help", argv[i])) {
			show_options(stdout);
			exit(0);
		}
		if (keyword("verbose", argv[i])) {
			verbose_flag = 1;
			continue;
		}
		if ((gotit=find("rms",argv[i])) != NULL) {
			rms=atof(gotit);
			continue;
		}
		if ((gotit=find("sky",argv[i])) != NULL) {
			sky=atof(gotit);
			continue;
		}
		if ((gotit=find("nrow",argv[i])) != NULL) {
			nrow=(int)(evaluate(gotit) + 0.5);
			continue;
		}
		if ((gotit=find("ncol",argv[i])) != NULL) {
			ncol=(int)(evaluate(gotit)+0.5);
			continue;
		}
		if ((gotit=find("rowgrad",argv[i])) != NULL) {
			rowgrad=atof(gotit);
			continue;
		}
		if ((gotit=find("colgrad",argv[i])) != NULL) {
			colgrad=atof(gotit);
			continue;
		}
		if ((gotit=find("gain",argv[i])) != NULL) {
			gain=atof(gotit);
			continue;
		}
		if ((gotit=find("seed",argv[i])) != NULL) {
			srand(atoi(gotit));
			continue;
		}
		if ((gotit=find("outfile",argv[i])) != NULL) {
			strcpy(fname,gotit);
			continue;
		}
		if ((gotit=find("star",argv[i])) != NULL) {
			nfield=sscanf(gotit,"%lf,%lf,%lf,%lf", 
				&s[nstar].height,
				&s[nstar].row,
				&s[nstar].col,
				&s[nstar].fwhm);
			if (NWIDTHS*s[nstar].fwhm > maxdist) {
				maxdist = NWIDTHS*s[nstar].fwhm;
			}
			if (nfield != 4) {
				error(1,"not enough fields in star (type help)");
			}
			nstar++;
			if (nstar == MAXSTAR) {
				error(1,"too many stars specified");
			}
			continue;
		}
		if ((gotit=find("field",argv[i])) != NULL) {
			nfield=sscanf(gotit,"%d,%lf,%lf,%lf", 
				&f.nstar,
				&f.max,
				&f.min,
				&f.fwhm);
			if (nfield != 4) {
				error(1,"not enough fields in field (type help)");
			}
			got_field=1;
			maxdist = NWIDTHS*f.fwhm;
			continue;
		}
		if ((gotit=find("grid",argv[i])) != NULL) {
			nfield=sscanf(gotit,"%d,%lf,%lf,%lf", 
				&f.nstar,
				&f.max,
				&f.min,
				&f.fwhm);
			if (nfield != 4) {
				error(1,"not enough fields in field (type help)");
			}
			got_grid = 1;
			maxdist = NWIDTHS*f.fwhm;
			continue;
		}
		if (keyword("photon",argv[i])) {
			photon = 1;
			continue;
		}
		fprintf(stderr,"unknown command '%s'\n",argv[i]);
		exit(1);
	}

	maxdist2 = maxdist*maxdist;
	negmaxdist = 0.0 - maxdist;

	if(got_field && got_grid){
		fprintf(stderr,"pstar: you can't make fields and grids together\n");
		exit(1);
	}
	if(got_field)
		nstar=make_field(s,&f,nrow,ncol);
	if(got_grid)
		nstar=make_grid(s,&f,nrow,ncol);
	

	fh=fits_open(fname, "w", &nrow, &ncol);

	fits_put_symbol(fh,"OBJECT","'artificial star'");
	sprintf(buf,"%.3f",sky);		fits_put_symbol(fh,"SKY",buf);
	sprintf(buf,"%.3f",rms);		fits_put_symbol(fh,"RMS",buf);
	sprintf(buf,"%.3f",gain);		fits_put_symbol(fh,"GAIN",buf);
	sprintf(buf,"%.3f",colgrad);	fits_put_symbol(fh,"COLGRAD",buf);
	sprintf(buf,"%.3f",rowgrad);	fits_put_symbol(fh,"ROWGRAD",buf);

	fits_put_symbol(fh,"PHOTON",photon?"T":"F");
	for(j=0; j<nstar; j++){
		sprintf(buf,"STAR%d",j);
		sprintf(buf2,"'%6.2f %6.2f %6.2f %6.2f'",
				s[j].height,
				s[j].row,
				s[j].col,
				s[j].fwhm);
		fits_put_symbol(fh,buf,buf2);
	}

	for(i=0; i<nrow; i++){
		if ((verbose_flag) && (i % 20 == 0)) {
			printf(" up to row %4d of %4d\r", i, nrow);
			fflush(stdout);
		}
		start=sky+rowgrad* (double)i;
		for(j=0; j<ncol; j++){
			val=start+colgrad*(double)j;
			for(k=0; k<nstar; k++){
				/* we assume all stars have same FWHM */
				dx = (double) i - s[k].row;
				dy = (double) j - s[k].col;
				if ((dx > negmaxdist) && (dx < maxdist) &&
				    (dy > negmaxdist) && (dy < maxdist)) {
					if (dx*dx + dy*dy < maxdist2) {
						val += s[k].height*dgauss(dx, dy,
								s[k].fwhm);
					}
				}
			}
			if(photon)
				val=poidev(gain*val,&idum)/gain;
			if(rms>0.)
				val+=(rms*gasdev(&idum));
			if(!warned && (val > (double)MAX_DATA_VAL)){
				fprintf(stderr,"warning: data exceeds %d\n",MAX_DATA_VAL);
				warned=1;
			}
			data[j]=(int16)floor(val+0.5);
		}
		fits_put_data(fh,i,0,data,ncol);
	}
	fits_close(fh);
	exit(0);
}

#define IGNORE	8.0			/* maximum distance to calculate */
#define FW_FACTOR 2.35482004503094938200

static double dgauss(dx,dy,s)
double dx,dy,s;
{
	double f;
	
	s/=FW_FACTOR;
	f=(dx*dx+dy*dy)/(2.0*s*s);
	if(f>IGNORE)
		return 0.;
	else
		return exp(-1.0*f);
}


static void 
show_options(fp)
FILE *fp;
{
	printf("pstar - generates artificial star image \n");
	printf("   [sky=]       - bias value  (default: 0.)\n");
	printf("   [rms=]       - additive r.m.s. noise (default: 0.)\n");
	printf("   [nrow=]      - number of rows in image (default:%d)\n",
				DEFAULT_SIZE);
	printf("   [ncol=]      - number of columns in image (default:%d) \n",
				DEFAULT_SIZE);
	printf("   [colgrad=]   - gradient per column (default: 0.)\n");
	printf("   [rowgrad=]   - gradient per row (default: 0.)\n");
	printf("   [photon]     - turns on photon noise (default: off)\n");
	printf("   [gain=]      - scale factor adus/photon (default: 1.0)\n");
	printf("   [outfile=]   - filename (default: %s%s)\n",
				DEFAULT_NAME,FITS_EXTENSION);
	printf("   [seed=]      - change random number seed value (default: 1)\n");
	printf("   [star=counts,row,col,fwhm]    - artificial star \n");
	printf("   [field=nstar,max,min,fwhm]    - field of artificial stars \n");
	printf("   [grid=nstar,max,min,fwhm]     - grid of artificial stars \n");
	printf("   [verbose]    - print progress through calculations \n");
}

/*	returns uniform random noise with zero mean and
	unit variance 
		Numerical Recipes Chapter 7
*/


static double gasdev(idum)
int *idum;
{
	static int iset;
	static double gset;
	double fac,r,v1,v2,val;

	if(iset == 0){
		do {
			v1=2.0*ran1(idum)-1.0;
			v2=2.0*ran1(idum)-1.0;
			r=v1*v1+v2*v2;
		} while( r>=1.0);
		fac=sqrt(-2.0*log(r)/r);
		gset=v1*fac;
		iset=1;
		val=v2*fac;
	}else {
		iset=0;
		val=gset;
	}
	return val;
}

/*	returns a Poisson deviate from the mean value of xm
	Stolen from Numerical Recipes
*/

static double poidev(xm, idum)
double xm;
int *idum;
{
	static double sq,alxm,g,oldm=(-1.0);
	double em,t,y;

	if (xm < 12.0){			/* use the direct method */
		if( xm !=oldm ) {
			oldm=xm;
			g=exp(-xm);
		}
		em =-1.0;
		t=1.0;
		do {			
			em+=1.0;
			t*= ran1(idum);
		} while(t>g);
	}else{
		if( xm!= oldm){
			oldm=xm;
			sq=sqrt(2.0*xm);
			alxm=log(xm);
			g=xm*alxm-gammln(xm+1.0);
		}
		do {
			do{
				y=tan(M_PI*ran1(idum));
				em=sq*y+xm;
			} while (em < 0.);
			em=floor(em);
			t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
		}while(ran1(idum)>t);
	}
	return em;
}

/*	Gamma function from Numerical Recipes */

static double gammln(xx)
double xx;
{
	double x,tmp,ser;
	static double cof[]={76.18009173,-86.50532033,24.0109822,
			-1.231739516, 0.120858003e-2,-0.536382e-5};
	int j;


	x=xx-1.0;
	tmp=x+5.5;
	tmp -= (x+0.5)*log(tmp);
	ser=1.0;
	for(j=0; j<=5; j++){
		x +=1.0;
		ser +=cof[j]/x;
	}
	return -tmp+ log(2.50662827*ser);
}

/*	poor uniform deviate between 0. and 1.0
*/

static double ran1(idum)
int *idum;
{
	return (double)rand()/(RAND_MAX);
}

/*	generates an array of f.nstar stars with peak values
	between f->max and f->min randomly picked in magnitude space

	the targets are sorted in row order and printed out to look
	like a file from 'stars'
*/


static int make_field(s,f,nrow,ncol)
struct star_type *s;
struct field_type *f;
int nrow,ncol;
{
	int i,idum;
	double max,min;

	if(f->max < f->min)
		error(1,"pstar: field maximum must be greater than minimum");
	if(f->min <= 0)
		error(1,"pstar: field minimum value must be greater than 0.");
	if(f->nstar >= MAXSTAR) 
		error(1,"pstar: field has too many stars");

	max=log(f->max);
	min=log(f->min);


	for(i=0; i<f->nstar; i++){
		s[i].row=ran1(&idum)*nrow;
		s[i].col=ran1(&idum)*ncol;
		s[i].fwhm=f->fwhm;
		s[i].height=exp((max-min)*ran1(&idum)+min);
	}

	qsort((void *) s,f->nstar,sizeof(struct star_type), (PFI) sort_by_row);

	for(i=0; i<f->nstar; i++)
		printf("%5d %7.2f   %7.2f   %5.1f   %6.3f   %6.3f   %6.3f \n",
					i,
					s[i].row,
					s[i].col,
					s[i].height,
					s[i].fwhm,
					0.,
					0.5);

	return f->nstar;
}


/*	generate a grid of stars with the brightest at lower rows
	and with equal magnitude steps */

static int make_grid(s,f,nrow,ncol)
struct star_type *s;
struct field_type *f;
int nrow,ncol;
{
	int i;
	double step,row,col,max,min;

	if(f->max < f->min)
		error(1,"pstar: grid maximum must be greater than minimum");
	if(f->min <= 0)
		error(1,"pstar: grid minimum value must be greater than 0.");
	if(f->nstar >= MAXSTAR) 
		error(1,"pstar: grid has too many stars");

	step=sqrt(nrow*ncol/(double)f->nstar+2.)	;	/* approximation */

	row=step/2;
	col=step/2;

	max=log(f->max);
	min=log(f->min);

	for(i=0; i<f->nstar; i++){
		s[i].row=row;
		s[i].col=col;
		s[i].fwhm=f->fwhm;
		s[i].height=exp((min-max)*((double)i/f->nstar)+max);
		col += step;
		if(col > (double) ncol){
			col=step/2.;
			row += step;
		}
	}

	for(i=0; i<f->nstar; i++)
		printf("%5d %7.2f   %7.2f   %5.1f   %6.3f   %6.3f   %6.3f \n",
					i,
					s[i].row,
					s[i].col,
					s[i].height,
					s[i].fwhm,
					0.,
					0.5);

	return f->nstar;
}

int sort_by_row(a,b)
struct star_type *a,*b;
{
	if( a->row > b->row)
		return 1;
	if( a->row < b->row)
		return -1;
	return 0;
}