#include "astrolib.h"
#include "fitting.h"

#ifdef WIN
#pragma warning (disable: 4018) // disable warning for comparing signed and unsigned
#pragma warning (disable: 4065) // disable warning for switch with default without case
#endif /* WIN */

/*
; create a nice string of celestial coordinates.
; ra and dec should be given in radians.
*/
char *coordstring(double ra, double dec, char *buffer) {
	static double rad_to_deg = 180.0 / M_PI;
	char sign = dec < 0 ? '-' : '+';
	double adec = fabs(dec);
	double rad = ra * rad_to_deg / 15.0;
	int rah = (int)rad;
	rad = (rad-rah) * 60.0;
	int ram = (int)rad;
	rad = (rad-ram) * 60.0;
	double de = adec * rad_to_deg;
	int deh = (int)de;
	de = (de-deh) * 60.0;
	int dem = (int)de;
	de = (de-dem) * 60.0;

	sprintf(buffer, "RA = %2ih %2im %5.2fs, DEC = %c%2id %2i' %4.1f\"", rah, ram, rad, sign, deh, dem, de);
	return buffer;
}

double ten(const double dd, const double mm, const double ss) {
	double result;

	result = fabs(dd) + fabs(mm) / 60.0 + fabs(ss) / 3600;
	if ((dd < 0.0) || (mm < 0.0) || (ss < 0.0)) result *= -1.0;

	return result;
}

void premat(const double equinox1, const double equinox2, double *r) {
	static double deg_to_rad = M_PI / 180.0;
	static double sec_to_rad = deg_to_rad / 3600;
	double T = 0.001*( equinox2 - equinox1);
	double ST = 0.001*( equinox1 - 2000.0);

//  Compute 3 rotation angles
	double A = sec_to_rad * T * (23062.181 + ST*(139.656 +0.0139*ST) + T*(30.188 - 0.344*ST+17.998*T));
	double B = sec_to_rad * T * T * (79.280 + 0.410*ST + 0.205*T) + A;
	double C = sec_to_rad * T * (20043.109 - ST*(85.33 + 0.217*ST) + T*(-42.665 - 0.217*ST -41.833*T));
	double sina = sin(A); double sinb = sin(B); double sinc = sin(C);
	double cosa = cos(A); double cosb = cos(B); double cosc = cos(C);

	r[0] =  cosa*cosb*cosc-sina*sinb;
	r[1] =  sina*cosb+cosa*sinb*cosc;
	r[2] =  cosa*sinc;
	r[3] = -cosa*sinb-sina*cosb*cosc;
	r[4] =  cosa*cosb-sina*sinb*cosc;
	r[5] = -sina*sinc;
	r[6] = -cosb*sinc;
	r[7] = -sinb*sinc;
	r[8] =  cosc;
}

void precess(double *ra, double *dec, const double equinox1, const double equinox2) {
	static double deg_to_rad = M_PI/180;

	double ra_rad = *ra * deg_to_rad;
	double dec_rad = *dec * deg_to_rad;

	double a = cos(dec_rad);
	double x[3];
  
	x[0] = a*cos(ra_rad);
	x[1] = a*sin(ra_rad);
	x[2] = sin(dec_rad);

// Use PREMAT function to get precession matrix from Equinox1 to Equinox2
	double r[9];
	premat(equinox1, equinox2, r);
	double x2[3];

	x2[0] = r[0] * x[0] + r[3] * x[1] + r[6] * x[2];
	x2[1] = r[1] * x[0] + r[4] * x[1] + r[7] * x[2];
	x2[2] = r[2] * x[0] + r[5] * x[1] + r[8] * x[2];

	ra_rad = atan2(x2[1],x2[0]);
	dec_rad = asin(x2[2]);
	*ra = ra_rad ; *dec = dec_rad;
	if (*ra < 0) *ra += 2.0*M_PI;

	char buffer[100];
	dp_output("Precessed coords are: %s\n", coordstring(*ra, *dec, buffer));
}

void precess(double rah, double ram, double ras, double ded, double dem, double des, const double equinox1, const double equinox2) {
	double ra = ten(rah, ram, ras) * 15.0;
	double dec = ten(ded, dem, des);

	precess(&ra, &dec, equinox1, equinox2);
}

bool find(Fits &image, Fits &result, double hmin, double fwhm, Fits &roundlim, Fits &sharplim, bool PRINT) {
	long i, _i;
	int maxbox = 13; // 	;Maximum size of convolution box in pixels 

	if (image.Naxis(0) != 2) {
		dp_output("find: ERROR - Image array (first parameter) must be 2 dimensional\n");
		return FALSE;
	}
	
	int n_x = image.Naxis(1);
	int n_y = image.Naxis(2);

	dp_output("Input Image Size is %i by %i\n", n_x, n_y);

	double radius = 0.637*fwhm; //;Radius is 1.5 sigma
	if (radius < 2.001) radius = 2.001;
	double radsq = radius * radius;
	int nhalf = int(radius);
	if (nhalf > (maxbox-1)/2) nhalf = (maxbox-1) / 2;
	int nbox = 2*nhalf + 1; //	;# of pixels in side of convolution box 
	int middle = nhalf; //          ;Index of central pixel

	int lastro = n_x - nhalf;
	int lastcl = n_y - nhalf;
	double sigsq = (fwhm / 2.35482) * (fwhm / 2.35482);
	Fits mask;
	mask.create(nbox, nbox, I1); //   ;Mask identifies valid pixels in convolution box 
	Fits c;
	c.create(nbox, nbox); //      ;c will contain Gaussian convolution kernel

	Fits dd, dd2, w;
	dd.create(nbox - 1, 1);
	for (_i = 0; _i < nbox-1; _i++) dd.r4data[_i] = _i + 0.5 - middle;
	dd2.copy(dd);
	dd2.mul(dd);
	w.copy(dd);
	for (_i = 0; _i < nbox-1; _i++) w.r4data[_i] = 1.0 - 0.5 * (fabs(w.r4data[_i]) - 0.5) / (middle - 0.5);
	int ir = (nhalf-1);
	if (ir < 1) ir = 1;

	Fits row2;
	row2.create(nbox, 1);
	for (_i = 0; _i < nbox; _i++) row2.r4data[_i] = (_i - nhalf) * (_i - nhalf);
	
	Fits temp;
	for (i = 0; i <= nhalf; i++) {
		temp.copy(row2);
		temp.add((float)i*(float)i);
		for (_i = 0; _i < nbox; _i++) {
			c.r4data[c.C_I(_i,nhalf-i)] = temp[_i];
			c.r4data[c.C_I(_i,nhalf+i)] = temp[_i];
		}
	}
	for (_i = 0; _i < c.Nelements(); _i++) {
		if (c.r4data[_i] <= radsq) mask.i1data[_i] = 1;
	}

	int pixels = 0;
	Fits good;

	good.create((int)mask.get_flux(), 1, I4);
	for (_i = 0; _i < c.Nelements(); _i++) {
		if (mask.i1data[_i] != 0) {
			good.i4data[pixels] = _i;
			pixels++;
		}
	}
	// Value of c are now equal to distance to center

	c *= mask;
	// Make c into a Gaussian kernel
	for (_i = 0; _i < pixels; _i++) c.r4data[good.i4data[_i]] = exp(-0.5 * c.r4data[good.i4data[_i]] / sigsq);
	double sumc = c.get_flux();
	Fits c1;
	c1.copy(c);
	c1.mul(c);
	double sumcsq = c1.get_flux() - sumc * sumc / pixels;
	sumc = sumc/pixels;
	for (_i = 0; _i < pixels; _i++) c.r4data[good.i4data[_i]] = (c.r4data[good.i4data[_i]] - sumc)/sumcsq;
	c1.copy(row2);
	for (_i = 0; _i < c1.Nelements(); _i++) c1.r4data[_i] = exp(-.5*c1.r4data[_i]/sigsq);
	double sumc1 = c1.get_flux()/nbox;
	Fits cc;
	cc.copy(c1);
	cc.mul(c1);
	double sumc1sq = cc.get_flux() - sumc1;
	for (_i = 0; _i < c1.Nelements(); _i++) c1.r4data[_i] = (c1.r4data[_i]-sumc1)/sumc1sq;
	sumc = w.get_flux();  //                       ;Needed for centroid computation

	double _f = 0.0;
	for (_i = 0; _i < pixels; _i++) _f += c.r4data[good.i4data[_i]] * c.r4data[good.i4data[_i]];
	dp_output("RELATIVE ERROR computed from FWHM2: %f\n", sqrt(_f));

	dp_output("Beginning convolution of image\n");

	Fits h;
	int _x, _y;
	h.copy(image);
	h.convol(c); //    ;Convolve image with kernel "c"
	for (_y = 0; _y < h.Naxis(2); _y++) {
		for (_x = 0; _x <= nhalf - 1; _x++) h.r4data[h.C_I(_x, _y)] = 0.0;
		for (_x = n_x - nhalf; _x <= n_x - 1; _x++) h.r4data[h.C_I(_x, _y)] = 0.0;
	}
	for (_x = 0; _x < h.Naxis(1); _x++) {
		for (_y = 0; _y <= nhalf - 1; _y++) h.r4data[h.C_I(_x, _y)] = 0.0;
		for (_y = n_y - nhalf; _y <= n_y - 1; _y++) h.r4data[h.C_I(_x, _y)] = 0.0;
	}

	dp_output("Finished convolution of image\n");

	mask.i1data[mask.C_I(middle,middle)] = 0; //From now on we exclude the central pixel
	pixels--; //so the number of valid pixels is reduced by 1
	good.create((int)mask.get_flux(), 1, I4);
	pixels = 0;
	for (_i = 0; _i < c.Nelements(); _i++) {
		if (mask.i1data[_i] != 0) {
			good.i4data[pixels] = _i;
			pixels++;
		}
	} //"good" identifies position of valid pixels
	Fits xx, yy, offset; // x and y coordinate of valid pixels relative to the center
	xx.copy(good);
	yy.copy(good);
	offset.copy(good);
	for (_i = 0; _i < xx.Nelements(); _i++) {
		xx.i4data[_i] = (good.i4data[_i] % nbox) - middle;
		yy.i4data[_i] = (int)(good.i4data[_i] / nbox) - middle;
		offset.i4data[_i] = yy.i4data[_i] * n_x + xx.i4data[_i];
	}

	int nfound;
	Fits index, index1, stars;
//SEARCH: 			    // Threshold dependent search begins here
    index.dpGE(h, hmin); // Valid image pixels are greater than hmin
    nfound = index.mask2index();
	for (_i = 0; _i < nfound; _i++) index.i4data[_i]--;
	if (nfound == 0) { // Any maxima found?
		dp_output("ERROR - No maxima exceed input threshold of %f\n", hmin);
		return FALSE;
	}

	Fits _h, __h;
	for (i = 0; i < pixels; i++) {
		_h.create(index.Nelements(), 1);
		__h.copy(_h);
        dp_output("%i %i\n", index.Nelements(), h.Nelements());
		for (_i = 0; _i < index.Nelements(); _i++) {
            _h.r4data[_i] = h.r4data[index.i4data[_i]];
			__h.r4data[_i] = h.r4data[index.i4data[_i] + offset.i4data[i]];
		}
        stars.dpGE(_h, __h);
        nfound = stars.mask2index();
		for (_i = 0; _i < nfound; _i++) stars.i4data[_i]--;
	    if (nfound <= 0) { // Do valid local maxima exist?
			dp_output("ERROR - No maxima exceed input threshold of %f\n", hmin);
			return FALSE;
		}
		index1.create(nfound, 1, I4);
		for (_i = 0; _i < nfound; _i++) index1.i4data[_i] = index.i4data[stars.i4data[_i]];
		index.copy(index1);
	}

	Fits ix, iy;
	ix.create(index.Nelements(), 1, I4);
	iy.copy(ix);

	for (_i = 0; _i < index.Nelements(); _i++) {
		ix.i4data[_i] = index.i4data[_i] % n_x; // X index of local maxima
		iy.i4data[_i] = index.i4data[_i] / n_x; // Y index of local maxima
	}
	int ngood = index.Nelements();
	dp_output("%i local maxima located above threshold\n", ngood);

	int nstar = 0, badround = 0, badsharp = 0, badcntrd = 0;
	result.create(5, ngood);
	Fits derivat, derivat2;
	double d, sharp1, dx, dy, around, sumd, sumxd, sumxsq, xcen, ycen;

//  Loop over star positions; compute statistics
	for (i = 0; i < ngood; i++) {
		temp.create(nhalf * 2 + 1, nhalf * 2 + 1);
		_i = 0;
		for (_x = ix[i] - nhalf; _x <= ix[i] + nhalf; _x++) {
			for (_y = iy[i] - nhalf; _y <= iy[i] + nhalf; _y++) {
				temp.r4data[_i] = image.ValueAt(image.C_I(_x, _y));
				_i++;
			}
		}
		d = h.r4data[h.C_I(ix[i],iy[i])];                  //"d" is actual pixel intensity        

//  Compute Sharpness statistic
		_h.copy(temp);
		_h.mul(mask);
		sharp1 = (temp.r4data[temp.C_I(middle,middle)] - (_h.get_flux())/pixels)/d;
		if ((sharp1 < sharplim[0]) || (sharp1 > sharplim[1])) {
			badsharp++;
			goto REJECT; //             ;Does not meet sharpness criteria
		}

//   Compute Roundness statistic
		dx = dy = 0.0;
		_h.collapse(temp, 1);
		_h.mul(c1);
		dx = _h.get_flux();

		_h.collapse(temp, 2);
		_h.mul(c1);
		dy = _h.get_flux();
		if ((dx <= 0) || (dy <= 0)) {
			badround++;
			goto REJECT; //           ;Cannot compute roundness
		}

		around = 2*(dx-dy) / ( dx + dy ); //    ;Roundness statistic
		if ((around < roundlim[0]) || (around > roundlim[1])) {
			badround++;
			goto REJECT; //           ;Does not meet roundness criteria
		}

// Find X centroid
		derivat.copy(temp);
		derivat.wrap(0,-1,0);
		derivat.sub(temp);
		derivat.extractRange(derivat2,middle-ir+1,middle+ir+1, 1,nbox-1, -1, -1);
		derivat.copy(derivat2);
		derivat2.collapse(derivat, 1);
		derivat.copy(derivat2);
		derivat.mul(w);
		sumd = derivat.get_flux();
		derivat.mul(dd);
		sumxd = derivat.get_flux();
		derivat.copy(w);
		derivat.mul(dd2);
		sumxsq = derivat.get_flux();

		if (sumxd >= 0.0) {
			badcntrd++;
			goto REJECT; //           ;Cannot compute X centroid
		}

		dx =sumxsq*sumd/(sumc*sumxd);
		if (fabs(dx) > nhalf) {
			badcntrd++;
			goto REJECT; //           ;X centroid too far from local X maxima
		}

		xcen = ix[i]-dx; //               ;Convert back to big image coordinates

// Find Y centroid                 
		derivat.copy(temp);
		derivat.wrap(-1,0,0);
		derivat.sub(temp);
		derivat.extractRange(derivat2, 1,nbox-1,middle-ir+1,middle+ir+1, -1, -1);
		derivat.copy(derivat2);
		derivat2.collapse(derivat, 2);
		derivat.copy(derivat2);
		derivat.mul(w);
		sumd = derivat.get_flux();
		derivat.mul(dd);
		sumxd = derivat.get_flux();
		derivat.copy(w);
		derivat.mul(dd2);
		sumxsq = derivat.get_flux();

		if (sumxd >= 0.0) {
			badcntrd++;
			goto REJECT; //           ;Cannot compute X centroid
		}

		dy =sumxsq*sumd/(sumc*sumxd);
		if (fabs(dy) > nhalf) {
			badcntrd++;
			goto REJECT; //           ;X centroid too far from local X maxima
		}

		ycen = iy[i]-dy; //               ;Convert back to big image coordinates

//  This star has met all selection criteria.  Print out and save results

		result.r4data[result.C_I(0, nstar)] = xcen + 1.0;
		result.r4data[result.C_I(1, nstar)] = ycen + 1.0;
		result.r4data[result.C_I(2, nstar)] = d;
		result.r4data[result.C_I(3, nstar)] = sharp1;
		result.r4data[result.C_I(4, nstar)] = around;

		if (PRINT) {
			if (nstar == 0) dp_output("   STAR      X        Y       FLUX      SHARP     ROUND\n");
			dp_output("%7i %8.1f %8.1f %9.1f %9.2f %9.2f\n", nstar+1, xcen+1.0, ycen+1.0, d, sharp1, around);
		}

		nstar++;
   
REJECT: ;
	}

	dp_output(" No. of sources rejected by SHARPNESS criteria: %i\n", badsharp);
	dp_output(" No. of sources rejected by ROUNDNESS criteria: %i\n", badround);
	dp_output(" No. of sources rejected by CENTROID  criteria: %i\n", badcntrd);

	if (nstar > 0) result.resize(5, nstar+1);

	return TRUE;
}

void abszissaGenerate(Fits &values, const Fits &data, int axis) {
	unsigned long i;
	double crval, crpix, cdelt;
	
	values.create(data.Naxis(axis), 1, R8);
	crval = data.getCRVAL(axis);
	crpix = data.getCRPIX(axis);
	cdelt = data.getCDELT(axis);

	for (i = 1; i <= data.Naxis(axis); i++) {
		values.r8data[i-1] = crval + ((double)i - crpix) * cdelt;
	}
}

bool fitcontinuum(Fits &cont, const Fits &x, const Fits &y) {
	Fits _err, cont2, fitter, newy, newx;		
	double deviation;
	unsigned int i, j;

	_err.create(x.Naxis(1), 1, R8);
	newy.copy(_err);
	newx.copy(_err);
	_err = 1.0;
	polyfit1d(fitter, x, y, _err, 1);
	cont.copy(x);
	cont *= fitter.ValueAt(1);
	cont += fitter.ValueAt(0);
	cont2.copy(y);
	cont2 -= cont;
	cont2.cblank();
	deviation = cont2.get_meddev();
	j = 0;
	for (i = 0; i < x.Nelements(); i++) {
		if (fabs(y.ValueAt(i) - cont.ValueAt(i)) < deviation) {
			j++;
			newy.setValue(y.ValueAt(i), j, 1);
			newx.setValue(x.ValueAt(i), j, 1);
		}
	}
	if (j < 10) {
		cont = 0.0;
		return FALSE;
	} else {
		newx.resize(j);
		newy.resize(j);
		_err.resize(j);
		polyfit1d(fitter, newx, newy, _err, 2);
		cont.copy(x);
		cont2.copy(x);
		cont2 *= cont2;
		cont2 *= fitter.ValueAt(2);
		cont *= fitter.ValueAt(1);
		cont += fitter.ValueAt(0);
		cont += cont2;
	}
	
	return TRUE;
}

/*
 cross-correlate an object with a template.
 The object and the template are rebinned to the same (smallest) dispersion
 and clipped to the overlapping spectral region.
 The variable option defines what should be done before
 doing the cross-correlation:
  option & 1: fit a continuum to the object and subtract
  option & 2: fit a continuum to the template and subtract
 Those values can be or'ed to perform both operations.
*/

bool fxcor(Fits &result, Fits &object, const Fits &templat, int option) {
	Fits xt, xo, yt, yo, x, cont;
	double xmin, xmax, xdelt;
//	unsigned long i;
	int objmin, objmax;

	xdelt = templat.getCDELT(1);
	if (object.getCDELT(1) > xdelt) xdelt = object.getCDELT(1);
	abszissaGenerate(xt, templat, 1);
	abszissaGenerate(xo, object, 1);

	xmin = xt.r8data[0];
	if (xo.r8data[0] > xmin) xmin = xo.r8data[0];
	xmax = xt.r8data[xt.Nelements() - 1];
	if (xo.r8data[xo.Nelements() - 1] < xmax) xmax = xo.r8data[xo.Nelements() - 1];

	yt.copy(templat);
	yt.norm();
	yt.rebin1d(xmin, xmax, xdelt);

	objmin = (int)(object.getCRPIX(1) + (xmin - object.getCRVAL(1)) / object.getCDELT(1) + 0.5);
	if (objmin < 1) objmin = 1;
	objmax = (int)(object.getCRPIX(1) + (xmax - object.getCRVAL(1)) / object.getCDELT(1) + 0.5);
	if (objmax > object.Naxis(1)) objmax = object.Naxis(1);

	object.extractRange(yo, objmin, objmax, 1, 1, 1, 1);
	yo.norm();
	yo.rebin1d(xmin, xmax, xdelt);

	abszissaGenerate(x, yo, 1);

	if (option & 1) { // continuum subtract object
		if (!fitcontinuum(cont, x, yo)) return FALSE;		
		yo -= cont;
	}
	if (option & 2) { // continuum subtract template
		if (!fitcontinuum(cont, x, yt)) return FALSE;		
		yt -= cont;
	}
	result.copy(yt);
	result.correl_real(yo);
	result.reass();

	return TRUE;
}

bool correlmap(Fits &result, Fits &cube, const Fits &templat, int method) {
	long x, y;
	int dispersionaxis = 3;
	double xmin, xmax, xdelt;
	Fits spectrum, correlation, xt, xo, yt, cont;
	int objmin, objmax;
	
	if (cube.Naxis(0) != 3) return FALSE;
	result.create(cube.Naxis(1), cube.Naxis(2), R4);
	
	dpProgress(-cube.Naxis(1), (const char *)("correlmap"));

// find out inclusive spectral range of object and template and minimum dispersion
	xdelt = templat.getCDELT(1);
	if (cube.getCDELT(dispersionaxis) > xdelt) xdelt = cube.getCDELT(dispersionaxis);
	abszissaGenerate(xt, templat, 1);
	abszissaGenerate(xo, cube, dispersionaxis);

	xmin = xt.r8data[0];
	if (xo.r8data[0] > xmin) xmin = xo.r8data[0];
	xmax = xt.r8data[xt.Nelements() - 1];
	if (xo.r8data[xo.Nelements() - 1] < xmax) xmax = xo.r8data[xo.Nelements() - 1];

// prepare template: normalize, rebin, and continuum subtract
	yt.copy(templat);
	yt.norm();
	yt.rebin1d(xmin, xmax, xdelt);
	abszissaGenerate(xt, yt, 1);
	if (!fitcontinuum(cont, xt, yt)) return FALSE;		
	yt -= cont;

// find out range for object
	objmin = (int)(cube.getCRPIX(dispersionaxis) + (xmin - cube.getCRVAL(dispersionaxis)) / cube.getCDELT(dispersionaxis) + 0.5);
	if (objmin < 1) objmin = 1;
	objmax = (int)(cube.getCRPIX(dispersionaxis) + (xmax - cube.getCRVAL(dispersionaxis)) / cube.getCDELT(dispersionaxis) + 0.5);
	if (objmax > cube.Naxis(dispersionaxis)) objmax = cube.Naxis(dispersionaxis);
	
	for (x = 1; x <= cube.Naxis(1); x++) {
		for (y = 1; y <= cube.Naxis(2); y++) {
			cube.extractRange(spectrum, x, x, y, y, objmin, objmax);
			spectrum.norm();
			spectrum.rebin1d(xmin, xmax, xdelt);
			if (!fitcontinuum(cont, xt, spectrum)) result.r4data[result.F_I(x, y)] = 0.0;
			else {	
				spectrum -= cont;
				correlation.copy(yt);
				correlation.correl_real(spectrum);
				correlation.reass();
				switch (method) {
					default:
						result.r4data[result.F_I(x, y)] = correlation.get_max();
						break;
				}
			}
		}
		if (FitsInterrupt) {
			dpProgress(0, "correlmap");
			return TRUE;
		}
		dpProgress(x, "correlmap");
	}
	
	return TRUE;
}

bool longslit(Fits &result, Fits &cube, int xcenter, int ycenter, double angle, int width, double opening_angle) {
    int newx, newy, sx, sy, z, y, i, owidth;
    double value, tan_opening = tan(opening_angle * M_PI / 180.);
	Fits slice, line, linesum;
	char key[255];
	
	if (cube.Naxis(3) < 2) return FALSE;
	newx = (int)(cube.Naxis(1) * 1.5);
	newy = (int)(cube.Naxis(2) * 1.5);
        if (newx > newy) newy = newx;
        if (newy > newx) newx = newy;

        sx = (int)((newx - cube.Naxis(1)) / 2);
        sy = (int)((newy - cube.Naxis(2)) / 2);
//	sx = (int)(cube.Naxis(1) * .25);
//	sy = (int)(cube.Naxis(2) * .25);
	xcenter += sx;
	ycenter += sy;
	
	if (!result.create(cube.Naxis(3), newy, R4)) return FALSE;
	result.CopyHeader(cube);
	result.SetFloatKey("CRVAL1", cube.getCRVAL(3));
	result.SetFloatKey("CRPIX1", cube.getCRPIX(3));
	result.SetFloatKey("CDELT1", cube.getCDELT(3));
	if (cube.GetStringKey("CTYPE3", key)) result.SetStringKey("CTYPE1", key);
	else result.DeleteKey("CTYPE1");
	if (cube.GetStringKey("CUNIT3", key)) result.SetStringKey("CUNIT1", key);
	else result.DeleteKey("CUNIT1");
	for (i = 2; i <= MAXNAXIS; i++) {
		sprintf(key, "CUNIT%i", i);
		result.DeleteKey(key);
		sprintf(key, "CTYPE%i", i);
		result.DeleteKey(key);
		sprintf(key, "CRPIX%i", i);
		result.DeleteKey(key);
		sprintf(key, "CRVAL%i", i);
		result.DeleteKey(key);
		sprintf(key, "CDELT%i", i);
		result.DeleteKey(key);
	}
	result.DeleteKey("CROTA1");
	result.DeleteKey("CROTA2");
	result.DeleteKey("CD1_1");
	result.DeleteKey("CD1_2");
	result.DeleteKey("CD2_1");
	result.DeleteKey("CD2_2");
	result.DeleteKey("CD3_3");
	sprintf(result.crtype, "");
	
	dpProgress(-cube.Naxis(3), "longslit");
	for (z = 1; z <= cube.Naxis(3); z++) {
		dpProgress(z, "longslit");
		cube.extractRange(slice, -1, -1, -1, -1, z, z);
		slice.resize(newx, newy);
		slice.ishift(sx, sy, 0);
		slice.rotate(angle, xcenter, ycenter, FALSE);
		slice.extractRange(linesum, xcenter, xcenter, -1, -1, -1, -1);

        for (y = 1; y <= newy; y++) {
            owidth = width + fabs(y - ycenter) * tan_opening * 2.; // TODO: Calculate opening
            if (owidth % 2 == 1) {
                for (i = 1; i <= (int)(owidth / 2); i++) {
                    value = slice.ValueAt(slice.at(xcenter-i, y));
                    linesum.setValue(linesum.ValueAt(linesum.at(y, 1)) + value, y, 1);
                    value = slice.ValueAt(slice.at(xcenter+i, y));
                    linesum.setValue(linesum.ValueAt(linesum.at(y, 1)) + value, y, 1);
//                    slice.extractRange(line, xcenter-i, xcenter-i, -1, -1, -1, -1);
//                    linesum += line;
//                    slice.extractRange(line, xcenter+i, xcenter+i, -1, -1, -1, -1);
//                    linesum += line;
                }
            } else {
                for (i = 1; i <= (int)(owidth / 2 - 1); i++) {
                    value = slice.ValueAt(slice.at(xcenter-i, y));
                    linesum.setValue(linesum.ValueAt(linesum.at(y, 1)) + value, y, 1);
                    value = slice.ValueAt(slice.at(xcenter+i, y));
                    linesum.setValue(linesum.ValueAt(linesum.at(y, 1)) + value, y, 1);
//                    slice.extractRange(line, xcenter-i, xcenter-i, -1, -1, -1, -1);
//                    linesum += line;
//                    slice.extractRange(line, xcenter+i, xcenter+i, -1, -1, -1, -1);
//                    linesum += line;
                }
                value = slice.ValueAt(slice.at(xcenter-(int)(owidth / 2), y));
                linesum.setValue(linesum.ValueAt(linesum.at(y, 1)) + value / 2., y, 1);
                value = slice.ValueAt(slice.at(xcenter+(int)(owidth / 2), y));
                linesum.setValue(linesum.ValueAt(linesum.at(y, 1)) + value / 2., y, 1);
//                slice.extractRange(line, xcenter-(int)(width / 2), xcenter-(int)(width / 2), -1, -1, -1, -1);
//                line /= 2.;
//                linesum += line;
//                slice.extractRange(line, xcenter+(int)(width / 2), xcenter+(int)(width / 2), -1, -1, -1, -1);
//                line /= 2.;
//                linesum += line;
            }
        }
        result.setRange(linesum, z, z, -1, -1, -1, -1);
    }
	
	return TRUE;
}

bool twodcut(Fits &result, Fits &image, int xcenter, int ycenter, double angle, int width) {
	int newx, newy, sx, sy, /*z,*/ i;
	Fits slice, line, linesum;
	char key[255];
	
        if (image.Naxis(2) < 2) return FALSE;
// handle case 0 or 180 degrees (cut along y-axis)
        if (angle == 0. || angle == 180.) {
            image.extractRange(line, xcenter - (int)(width / 2), xcenter + (int)(width / 2), -1, -1, -1, -1);
            line.setType(R8);
            if (width % 2 == 0) {
                for (i = 0; i < line.Naxis(2); i++) {
                    line.r8data[line.C_I(0, i)] /= 2.;
                    line.r8data[line.C_I(width, i)] /= 2.;
                }
            }
            if (width < 2) {
                result.copy(line);
            } else {
                result.collapse(line, 1);
            }
            if (angle == 180.) result.flip(1);

            return TRUE;
        }
// handle case 90 or 270 degrees (cut along x-axis)
        if (angle == 90. || angle == 270.) {
            image.extractRange(line, -1, -1, ycenter - (int)(width / 2), ycenter + (int)(width / 2), -1, -1);
            line.setType(R8);
            if (width % 2 == 0) {
                for (i = 0; i < line.Naxis(1); i++) {
                    line.r8data[line.C_I(i, 0)] /= 2.;
                    line.r8data[line.C_I(i, width)] /= 2.;
                }
            }
            result.collapse(line, 2);
            if (angle == 270.) result.flip(1);

            return TRUE;
        }

	newx = (int)(image.Naxis(1) * 1.5);
	newy = (int)(image.Naxis(2) * 1.5);
        if (newx > newy) newy = newx;
        if (newy > newx) newx = newy;

        sx = (int)((newx - image.Naxis(1)) / 2);
        sy = (int)((newy - image.Naxis(2)) / 2);
	xcenter += sx;
	ycenter += sy;
	
	result.CopyHeader(image);
	result.SetFloatKey("CRVAL1", 1.0);
	result.SetFloatKey("CRPIX1", 1.0);
	result.SetFloatKey("CDELT1", (image.getCDELT(1) + image.getCDELT(2)) / 2.0);
	for (i = 2; i <= MAXNAXIS; i++) {
		sprintf(key, "CUNIT%i", i);
		result.DeleteKey(key);
		sprintf(key, "CTYPE%i", i);
		result.DeleteKey(key);
		sprintf(key, "CRPIX%i", i);
		result.DeleteKey(key);
		sprintf(key, "CRVAL%i", i);
		result.DeleteKey(key);
		sprintf(key, "CDELT%i", i);
		result.DeleteKey(key);
	}
	result.DeleteKey("CROTA1");
	result.DeleteKey("CROTA2");
	result.DeleteKey("CD1_1");
	result.DeleteKey("CD1_2");
	result.DeleteKey("CD2_1");
	result.DeleteKey("CD2_2");
	sprintf(result.crtype, "");
	
	slice.copy(image);
	slice.resize(newx, newy);
	slice.ishift(sx, sy, 0);
	slice.rotate(angle, xcenter, ycenter, FALSE);
	slice.extractRange(linesum, xcenter, xcenter, -1, -1, -1, -1);
	if (width % 2 == 1) {
		for (i = 1; i <= (int)(width / 2); i++) {
			slice.extractRange(line, xcenter-i, xcenter-i, -1, -1, -1, -1);
			linesum += line;
			slice.extractRange(line, xcenter+i, xcenter+i, -1, -1, -1, -1);
			linesum += line;
		}
	} else {
		for (i = 1; i <= (int)(width / 2 - 1); i++) {
			slice.extractRange(line, xcenter-i, xcenter-i, -1, -1, -1, -1);
			linesum += line;
			slice.extractRange(line, xcenter+i, xcenter+i, -1, -1, -1, -1);
			linesum += line;
		}
		slice.extractRange(line, xcenter-(int)(width / 2), xcenter-(int)(width / 2), -1, -1, -1, -1);
		line /= 2.;
		linesum += line;
		slice.extractRange(line, xcenter+(int)(width / 2), xcenter+(int)(width / 2), -1, -1, -1, -1);
		line /= 2.;
		linesum += line;
	}
	if (!result.copy(linesum)) return FALSE;
	result.CopyHeader(image);
	result.SetFloatKey("CRVAL1", 1.0);
	result.SetFloatKey("CRPIX1", 1.0);
	result.SetFloatKey("CDELT1", (fabs(image.getCDELT(1)) + fabs(image.getCDELT(2))) / 2.0);
	if (result.getCDELT(1) == 0.0) result.SetFloatKey("CDELT1", 1.0);
	for (i = 2; i <= MAXNAXIS; i++) {
		sprintf(key, "CUNIT%i", i);
		result.DeleteKey(key);
		sprintf(key, "CTYPE%i", i);
		result.DeleteKey(key);
		sprintf(key, "CRPIX%i", i);
		result.DeleteKey(key);
		sprintf(key, "CRVAL%i", i);
		result.DeleteKey(key);
		sprintf(key, "CDELT%i", i);
		result.DeleteKey(key);
	}
	result.DeleteKey("CROTA1");
	result.DeleteKey("CROTA2");
	sprintf(result.crtype, "");
	
	
	return TRUE;
}

Fits &primes(int k, Fits &rv) {
	rv.create(k, 1, I4);
	
	if (k == 1) {
		rv.i4data[0] = 2;
		return rv;
	}
	
	rv.i4data[0] = 2;
	unsigned long n = 3;
	unsigned long count = 1;
	rv.i4data[count] = 3;

	count++;

	while(count < k) {
		n = n + 2;

		for (int ip = 1; ip < count; ip++) {
			double q = n / rv.i4data[ip];
			int r = n % rv.i4data[ip];
			if (r == 0) {
				ip = count + 1;    // n is not prime.
			} else {
				if (q <= rv.i4data[ip]) {  // n is prime.
					rv.i4data[count] = n;
					count++;               // compute next prime.
					ip = count + 1;
				}
			}
		}
	}
	return rv;
}