//-----------------------------------------------------------------------------
// flatmapping.c  V7.1  S.F.Fulghum  1 September 2005
//
// ............................................................................
// Notes from 10 November 2007 when I finally posted the code
//
// email: steve@stevefulghum.com
// but I have probably forgotten how that part of the code works by now
// This code is something like dinosaur DNA
// you may know every letter... but without a dinosaur to lay the egg
// you can't build one.
// Mr. Koi himself is my own image taken in Kyoto.
// The galaxies are public domain Hubble Space Telescope image files
// which were rotated and stretched to first make them appear circular.
// All of the images were converted to the Photoshop RAW format
// of uncompressed binary numbers to make the I/O easier in C.
// The output is also in RAW format but was converted to TIF
// in Photoshop before printing.
//
// The original image was to have the right perspective if it was
// 15x10 inches in size, held by hand and normal to the line of sight
// from one's galactic eye looking down Mr. Koi
// floating in his galactic pond at one's galactic feet.
// Later the image was scaled larger so as not to waste
// the magnificent detail in the Hubble images of the galaxies.
// 15 inches at 360 DPI results in the original 5400 pixel width.
// Perspective lost out to detail as the image was scaled up
// to 6480 at 240 DPI. 
//
// this program compiles, but some variable are unreferenced
// because they were used in error checking or plotting grid lines
// and have been commented out 
// ............................................................................
// 
// 6480x4320@240 with perspective of 5400x3600@240
//
// adds saturation to the addition
// 6.0 adds black hole to distortion
// V6.1 is exponential hole
// V7.0 writes to more pixels (M51 doubled to minimize mapping error)
//      5400x3600=15x10@360DPI=22.5x15@240DPI
//		6480x4320=27x18@240DPI
//      
// adds galaxies to input image
// from mapping V12.0 25 April
// this modification starts with an image on the flat surface of the pool
// The idea is to modfify galaxy images until they are circular
// (easy since they are seen at a distance without perspective)
// rotate them appropriately in Photoshop
// clip them to a square format... place that square format on the flat surface
// add the ripple function
// then map them onto the picture plane where they will be foreshortened ellipes
// 
// The trick will be keeping the intensity of the mapping uniform between
// different pixel sizes of galaxies. High density images will appear brighter
// if some weighting is not applied to the individual images being mapped.
//
// The intent is to simulate waves in a thin field of stars.
// From an angled viewpoint the stars should appear more dense in some areas
// than in others, depending on the angle of intersection of the line of sight
// to the wave surface. 
// The final mapping will involve a sinusoidal surface containing the
// original pixels (by adding a patterned z offset to the 2-D image).
//
// Mapping proceeds by mapping each corner of a pixel on the flat surface
// onto the final image plane where it will be a rectangle in general
// which either fits within one pixel (if the density of pixels in the flat
// image is very high) or more generally among four or more pixels in the
// final image. The intensity of the flat pixel is subdivided appropriately.
//
// The final image should print at 360 pixels per inch onto a 10X15 inch region
// on the final image. This leaves 1.5 inches on all edges to keep the final
// image flat when matted. This means we need 3600x5400 final pixels (58 MB).
//
// Each galaxy to be mapped is treated sequentially with the same function
// This will allow as many images to be mapped as necessary by just copying
// the function with new inputs. Memory space will be freed between each.
//
// geometry: eye position is at x=0,y=0, z=H in flat space
//           x is left right, y increases with distance along flat space
//           picture is at distance D from eye, perpendicular to line of sight
//           to its center, which is off of the vertical by angle Td
//
// A function "flatmap" takes as input variables
// 1: a filename for a picture to be mapped
// 2: an intensity scale factor
// 3: the upper left corner x in flat space
// 4: the upper left corner y in flat space
// 5: the image size in pixels, ipx
// 6: the image size in pixels, ipy
// 7: the square pixel dimension in flat space, delta
//
//
//
// Global Variables

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

//old 5400x3600 values
//#define XPM1 5599
//#define YPM1 3799
//#define XPIX 5400
//#define YPIX 3600
//#define DPI 360
//#define OFF 100


// hopefully 432DPI will scale all the proportions correctly
// the intensity of each galaxy may need to be adjusted
// since it will be divided over more pixels

#define XPM1 6690
#define YPM1 4519
#define XPIX 6480
#define YPIX 4320
#define DPI 432
#define OFF 100

struct PIC
{
	double x;	// 2D position on picture measured from center
	double y;	// 2D position
	double ix;	// float integer in pixel space
	double iy;	// float integer
};

double **dmatrix(long nrl, long nrh, long ncl, long nch);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
void nrerror(char error_text[]);
void flatmap(char *flatimage, double Iscale, double fulx, double fuly,
			 int ifpx, int ifpy, double delta);
struct PIC map(double sx, double sy);

//=============================================================================
// Global Variables for projection

	double sint,cost;   // sin and cos of angle of view (fixed)
	double ex,ey,ez;    // position of eye in 3D space
	double cx,cy,cz;    // center of picture in 3-D space
	double Rcenx,Rceny; // 3D center of ripples on surface at z=0
	double Rwave,Ramp;  // wavelength and amplitude of ripple
	double R0;          // radius at which wave begins to oscillate
	double Rx,Ry;       // half widths of picture in 2D space
	double Step;        // size of pixel in 2D space
	double **outR;      // the red channel (double pointer to pixel)
	double **outB;      // the blue channel
	double **outG;      // the green channel
	double Isat;        // saturation parameters for summing
	double B,Bmin;      // horizontal black hole depth
	double xb,yb;       // black hole position
	double rb0,Rb;      // black hole inner and outer radius  
	double Sx,Sy,Sb;    // streamline strength into black hole

//*****************************************************************************
void main(void)
{

	double conpi,con2pi;  // useful constants
	double xi,yi;         // initial 2D position on image from center
	double bx,by,bz;      // initial position on pic in 3D space
	double xf,yf;         // 2D remapped x,y
	double sx,sy;         // 2D projected position on surface
	double sz;            // 3D z offset due to waves above xs,ys
	double fx,fy,fz;      // 3D mapped 3D intersection
	double t;             // angle in radians
	double fixm,fiym;
	double dtemp,Rtemp,Gtemp,Btemp;
	double Lbase,Lgalaxy;
	double xpixels,ypixels;
	

	double D;    // distance to image
	double H;    // height of eye above surface
	double Td;   // angle of obsrvation from vertical in degrees

	long satcount;
	short issat;

	int ix,iy;
	int igx,igy;
	int ixm,iym;
	double dgx,dgy;

	unsigned char ucval;
	unsigned short usval;

	FILE *fpbase;
	FILE *fpout;

	struct PIC grid;

	conpi=3.141592653589793238;
	con2pi=2.0*conpi;

// input parameters
// one good fit is:
// D=16  H=60  Td=38.5  Rwave=4.8
// xi=-0.1297 yi=0.825
//
// the 10x15 image requires D=28 instead of 16 to keep the perpsective
// the same as in the smaller, handheld image

	D= 28.0;
	H= 60.0;
	Td=38.5;

	B=-12.0;
	Bmin=-4.5;
	xb=-1.5;
	yb=47.0;
	rb0=0.2;
	Rb=4.0;
	Sb=1.0;
	Sx=1.5;
	Sy=1.4; // sets strength of streamline distortion

// b parameters: Isat=65000, needs to be lower to raise spiral arm visibility
// c parameters: Isat=50000
// d parameters: Isat=40000

 	Isat=40000.0;

//-----------------------------------------------------------------------------
// input or calculate global variable for rest of program

	Rwave=3.5;
	Ramp=0.5;
	R0=3.0*Rwave;   // radius at which wave begins to oscillate
	ex=0.0;
	ey=0.0;
	ez=H;
	t=conpi*Td/180.0;
	sint=sin(t);
	cost=sqrt(1.0-sint*sint);
	cx=0.0;
	cy=D*sint;
	cz=H-D*cost;

	printf("3D picture center: cx=%8.5f  cy=%8.5f  cz=%8.5f\n",cx,cy,cz);

//-----------------------------------------------------------------------------

	xpixels=(double)XPM1 +1.0;
	ypixels=(double)YPM1 +1.0;
	Step=1.0/(double)DPI;

	Rx=0.5*xpixels*Step;
	Ry=0.5*ypixels*Step;


// position of picture center from input

	sy=H*cy/(H-cz);
	sx=0.0;

	printf("\npicture center surface projection sy=%9.6f\n",sy); 
	printf("picture center surface projection sx=%9.6f\n",sx);

// calculate a center for ripples based on ellipse map
// center of ripples at picture position x=-0.1279, y=+0.7281
	
	xi= 0.30;
	yi= 0.10;

	bx=xi;
	by=cy+yi*cost;
	bz=cz+yi*sint;
	sy=H*by/(H-bz);
	sx=bx*sy/by;
	sz=0.0;

	Rcenx=sx;
	Rceny=sy;

	printf("\nripple center surface projection sy=%9.6f\n",sy);
	printf("ripple center surface projection sx=%9.6f\n",sx);

	fx=( sx*cy*cy+sx*(cz-ez)*(cz-ez) )/( cy*sy+(sz-ez)*(cz-ez) );                                          
	fy=( sy*cy*cy+sy*(cz-ez)*(cz-ez) )/( cy*sy+(sz-ez)*(cz-ez) );
	fz=( cz*(cz-ez)*(sz-ez)+cy*cy*(sz-ez)+cy*sy*ez )/( (cz-ez)*(sz-ez)+cy*sy );

	xf=fx;
	yf=(fy-cy)/cost;

	printf("\nripple center picture mapping xf=%9.6f\n",xf);
	printf("ripple center picture mapping yf=%9.6f\n",yf);

	fixm= ((xf+Rx)/Step) -0.5;
	fiym= ((Ry-yf)/Step) -0.5;
	
	printf("\nripple center pixel fixm=%9.6f\n",fixm);
	printf("ripple center pixel fiym=%9.6f\n",fiym);

//-----------------------------------------------------------------------------
// calculate pixel position of black hole

	sx=xb;
	sy=yb;
	sz=0.0;

	fx=( sx*cy*cy+sx*(cz-ez)*(cz-ez) )/( cy*sy+(sz-ez)*(cz-ez) );                                          
	fy=( sy*cy*cy+sy*(cz-ez)*(cz-ez) )/( cy*sy+(sz-ez)*(cz-ez) );
	fz=( cz*(cz-ez)*(sz-ez)+cy*cy*(sz-ez)+cy*sy*ez )/( (cz-ez)*(sz-ez)+cy*sy );

	xf=fx;
	yf=(fy-cy)/cost;

	printf("\nblack hole picture mapping xf=%9.6f\n",xf);
	printf("black hole picture mapping yf=%9.6f\n",yf);

	fixm= ((xf+Rx)/Step) -0.5;
	fiym= ((Ry-yf)/Step) -0.5;
	
	printf("\nblack hole pixel fixm=%9.6f\n",fixm);
	printf("black hole pixel fiym=%9.6f\n",fiym);


//-----------------------------------------------------------------------------
// input the image


	fpout=fopen("mrkoi_6480_4320_240_f.raw","wb");
	if (fpout==NULL) {
		printf("Cannot open output image file\n");
		exit(1);
	}

// small final image base file
//	fpbase=fopen("koi_base_5400_3600.raw","rb");

	fpbase=fopen("koi_base_6480_4320.raw","rb");
	if (fpbase==NULL) {
		printf("Cannot open input image file\n");
		exit(1);
	}
 
	printf("files found\n");
	printf("allocating space\n");

// zero the initial image storage space

	outR=dmatrix(0,XPM1,0,YPM1);
	printf("outR done\n");
	outG=dmatrix(0,XPM1,0,YPM1);
	printf("outG done\n");
	outB=dmatrix(0,XPM1,0,YPM1);
	printf("outB done\n");

	printf("space allocated\n");

	printf("zeroing %6d memory rows\n",YPM1);
	for(iy=0;iy<=YPM1;iy++) {
		printf("%5d\b\b\b\b\b",iy);
		for(ix=0;ix<=XPM1;ix++) {
			outR[ix][iy]=0.0;
			outG[ix][iy]=0.0;
			outB[ix][iy]=0.0;
		}
	}

// "b" parameters
// 3370 a bit dim, M51 and 1300 ok, need to lower saturation to bring up
// spiral arm structure, will need to raise intensity to compensate
// b Isat=65000
// flatmap("NGC1300_7364_6487.raw",60.0,-16.0,51.0,7364,6487,0.0034);
// flatmap("NGC3370_12296_12296.raw",160.0,-40.0,87.0,12296,12296,0.004);
// flatmap("M51_10688_9934.raw",70.0,6.5,64.0,10688,9934,0.0025);

// c parameters with Isat=50000
// flatmap("NGC1300_7364_6487.raw",66.0,-16.0,51.0,7364,6487,0.0034);
// flatmap("NGC3370_12296_12296.raw",200.0,-40.0,87.0,12296,12296,0.004);
// flatmap("M51_10688_9934.raw",77.0,6.5,64.0,10688,9934,0.0025);


// d parameters with Isat=40000
// flatmap("NGC1300_7364_6487.raw",75.0,-16.0,51.0,7364,6487,0.0034);
// flatmap("NGC3370_12296_12296.raw",230.0,-40.0,87.0,12296,12296,0.004);
// flatmap("M51_10688_9934.raw",85.0,6.5,64.0,10688,9934,0.0025);

// e parameters: an Isat=20000 emphasizes a Moire' pattern in the galaxies
// drop back to 40000 (d) but add 8% to NGC1300 for (e)
// flatmap("NGC1300_7364_6487.raw",85.0,-16.0,51.0,7364,6487,0.0034);
// flatmap("NGC3370_12296_12296.raw",230.0,-40.0,87.0,12296,12296,0.004);
// flatmap("M51_10688_9934.raw",85.0,6.5,64.0,10688,9934,0.0025);
// this value for NGC1300 brightens it without losing detail in galaxy center

// f parameter: Isat=40000, pick up brightness of M51 a bit, the best so far
// flatmap("NGC1300_7364_6487.raw",85.0,-16.0,51.0,7364,6487,0.0034);
// flatmap("NGC3370_12296_12296.raw",230.0,-40.0,87.0,12296,12296,0.004);
// flatmap("M51_10688_9934.raw",90.0,6.5,64.0,10688,9934,0.0025);


	printf("running flatmap\n");

 	flatmap("NGC1300_7364_6487.raw",85.0,-16.0,51.0,7364,6487,0.0034);

	flatmap("NGC3370_12296_12296.raw",230.0,-40.0,87.0,12296,12296,0.004);

	flatmap("M51_10688_9934.raw",90.0,6.5,64.0,10688,9934,0.0025);


	printf("flatmap complete, writing output image\n");

//=============================================================================
// find peak pixel value of galaxies and how many pixels are saturated
/*
	printf("finding peak pixel value\n");
	satcount=0;
	dtemp=0.0;
	for(iy=OFF;iy<YPIX+OFF;iy++) {
		printf("%4d of %4d\b\b\b\b\b\b\b\b\b\b\b\b",iy,YPM1);
		for(ix=OFF;ix<XPIX+OFF;ix++) {
			issat=0;
			if (dtemp<outR[ix][iy]) dtemp=outR[ix][iy];
			if (outR[ix][iy]>65535.0) issat=1;
			if (dtemp<outG[ix][iy]) dtemp=outG[ix][iy];
			if (outR[ix][iy]>65535.0) issat=1;
			if (dtemp<outB[ix][iy]) dtemp=outB[ix][iy];
			if (outR[ix][iy]>65535.0) issat=1;
			if (issat) satcount+=1;
		}
	}
	printf("float peak pixel value = %9.1f\n",dtemp);

	printf("saturated pixels = %8d\n",satcount);
*/
//=============================================================================
// read in base image (in 8-bit mode).. scale up to 16 bits afterwards)
// and add it to above result

 
	printf("reading base image\n");
	for(iy=OFF;iy<YPIX+OFF;iy++) {
		printf("%4d of %4d\b\b\b\b\b\b\b\b\b\b\b\b",iy,YPM1);
		for(ix=OFF;ix<XPIX+OFF;ix++) {
			fread(&ucval,sizeof(unsigned char),1,fpbase);
			Rtemp=(double)ucval*54000.0/255.0;
			fread(&ucval,sizeof(unsigned char),1,fpbase);
			Gtemp=(double)ucval*54000.0/255.0;
			fread(&ucval,sizeof(unsigned char),1,fpbase);
			Btemp=(double)ucval*54000.0/255.0;
			if (outR[ix][iy]<Rtemp) outR[ix][iy]=Rtemp;
			if (outG[ix][iy]<Gtemp) outG[ix][iy]=Gtemp;
			if (outB[ix][iy]<Btemp) outB[ix][iy]=Btemp;
		}
	}
	printf("\n");
 
// find peak pixel value of galaxies and how many piels are saturated
/*
	printf("finding peak pixel value after base is added\n");
	satcount=0;
	dtemp=0.0;
	for(iy=OFF;iy<YPIX+OFF;iy++) {
		printf("%4d of %4d\b\b\b\b\b\b\b\b\b\b\b\b",iy,YPM1);
		for(ix=OFF;ix<XPIX+OFF;ix++) {
			issat=0;
			if (dtemp<outR[ix][iy]) dtemp=outR[ix][iy];
			if (outR[ix][iy]>65535.0) issat=1;
			if (dtemp<outG[ix][iy]) dtemp=outG[ix][iy];
			if (outR[ix][iy]>65535.0) issat=1;
			if (dtemp<outB[ix][iy]) dtemp=outB[ix][iy];
			if (outR[ix][iy]>65535.0) issat=1;
			if (issat) satcount+=1;
		}
	}
	printf("float peak pixel value = %7.1f\n",dtemp);
	printf("saturated pixels = %8d\n",satcount);
*/
// write grid onto image for tests
/* 
	printf("write grids\n");

	for (dgy=20.0;dgy<=50.0;dgy+=0.5) {
		for (dgx=-10.0;dgx<=10.0;dgx+=0.01) {
			grid=map(dgx,dgy);
			ixm=(int)grid.ix;
			iym=(int)grid.iy;
			if ((ixm<3)||(iym<3)||(ixm>XPM1-3)||(iym>YPM1-3)) continue;
			outR[ixm][iym]=65500.0;
			outG[ixm][iym]=65500.0;
			outB[ixm][iym]=65500.0;
		}
	}
	for (dgx=-15.0;dgx<=15.0;dgx+=0.5) {
		for (dgy=20.0;dgy<=60.0;dgy+=0.01) {
			grid=map(dgx,dgy);
			ixm=(int)grid.ix;
			iym=(int)grid.iy;
			if ((ixm<3)||(iym<3)||(ixm>XPM1-3)||(iym>YPM1-3)) continue;
			outR[ixm][iym]=65500.0;
			outG[ixm][iym]=65500.0;
			outB[ixm][iym]=65500.0;
		}
	}
*/ 
// write modified image with interleaved RAW format

	printf("writing image to file\n");
	for(iy=OFF;iy<YPIX+OFF;iy++) {
		printf("%4d\b\b\b\b\b",iy);
		for(ix=OFF;ix<XPIX+OFF;ix++) {
			if(outR[ix][iy]>65500.0) outR[ix][iy]=65500.0;
			usval=(unsigned short)(outR[ix][iy]+0.000000001);
 			fwrite(&usval,sizeof(unsigned short),1,fpout);

			if(outG[ix][iy]>65500.0) outG[ix][iy]=65500.0;
			usval=(unsigned short)(outG[ix][iy]+0.000000001);
 			fwrite(&usval,sizeof(unsigned short),1,fpout);

			if(outB[ix][iy]>65500.0) outB[ix][iy]=65500.0;
			usval=(unsigned short)(outB[ix][iy]+0.000000001);
			fwrite(&usval,sizeof(unsigned short),1,fpout);
		}
	}
	fclose(fpout);

	printf("image written\n");

//-----------------------------------------------------------------------------
// free memory

	free_dmatrix(outR,0,XPM1,0,YPM1);
	free_dmatrix(outG,0,XPM1,0,YPM1);
	free_dmatrix(outB,0,XPM1,0,YPM1);

	printf("finis \n\n");


}
//*****************************************************************************
// functions
// map a flat surface point to a 2D picture point
// this includes calculating the vertical offset due to waves
// wave parameters are global
// black hole parameters are global

struct PIC map(double sx, double sy)
{
	double r;			// temps for wave calculation
	double rb;			// distance to black hole
	double sz;			// projected surface 3D position with offset
	double bsz;			// black hole vertical offset
	double bsx,bsy;		// horizontal offsets due to black hole
	double mx,my,mz;	// mapped 3D position on picture 
	struct PIC m;		// mapped 2D point for return value
	double sval;		// the streamline offset in the +y direction
	double sintheta,costheta;

// clip sy at zero to prevent blowups

	if (sy<0.0) sy=0.0;
	bsx=0.0;
	bsy=0.0;
	bsz=0.0;
	sz=0.0;


// add black hole horizontal offset and vertical (negative) offset

	rb=sqrt((sx-xb)*(sx-xb)+(sy-yb)*(sy-yb));
	costheta=(yb-sy)/rb;
	sintheta=(xb-sx)/rb;

	if (rb>rb0) bsz=B*exp(-(rb-rb0)/Rb);
	else bsz=B;
	if (bsz<Bmin) bsz=Bmin;

	sval=(rb/Sb)*(1.0-exp(-rb/Sb))*exp(-rb/(6.0*Sb));	
	bsx=Sx*sval*sintheta;	
	bsy=Sy*sval*costheta;

	sx+=bsx;
 	sy+=bsy;

// calculate ripple offset based on global Rwave,Ramp,Rcenx,Rceny
 
	r=sqrt((sx-Rcenx)*(sx-Rcenx) + (sy-Rceny)*(sy-Rceny));
	sz=Ramp*sin(2.0*3.14159265*r/Rwave);
	if (r<=R0) sz*=((r/R0)*(r/R0)*(r/R0)*(r/R0));
	else sz*=(R0/r);

	sz+=bsz;

// calculate 3D picture projection of offset ripple point

	mx=( sx*cy*cy+sx*(cz-ez)*(cz-ez) )/( cy*sy+(sz-ez)*(cz-ez) );                                          
	my=( sy*cy*cy+sy*(cz-ez)*(cz-ez) )/( cy*sy+(sz-ez)*(cz-ez) );
	mz=( cz*(cz-ez)*(sz-ez)+cy*cy*(sz-ez)+cy*sy*ez )/( (cz-ez)*(sz-ez)+cy*sy );

//	printf("3D mapped position mx=%9.6f  my=%9.6f  mz=%9.6f\n",mx,my,mz);
// note cost is a Global variable

	m.x=mx;
	m.y=(my-cy)/cost;

// convert to pixel space

	m.ix=(m.x + Rx)/Step;
	m.iy=(Ry - m.y)/Step;

//	printf("mapped pixel position: ix=%9.3f  iy=%9.3f\n",m.ix,m.iy);

	return m;

}

//=============================================================================
// map a flat image onto the final picture adding the surface ripple
// read the input image 1 pixel at a time so that storage is not necessary
// (we will try this... it may be more efficient to read it in all at once)


void flatmap(char *flatimage, double Iscale, double fulx,
			 double fuly, int ifx, int ify, double delta)
{
	struct PIC mul,mur,mll,mlr;	// mapped 2D corners of flat pixel on picture
	FILE *fpin;
	int ix,iy;			// pixel index of flat image
	double dix,diy;		//double precision ix,iy
	unsigned char ucval;// file value
	double R,G,B;		// RGB value of current pixel
	double I;			// current intensity of pixel
	double ulx,uly,urx,ury,llx,lly,lrx,lry;	// corners of current pixel
	int ixm,iym;		// integer values of upper left corner of mapped pixel
	int ox,oy,oxmax,oxmin,oymax,oymin,oxoy,oxoymax,oxoymin;
	int o1,o2,o3,o4,o5,o6,o7,o8,o9,o10,o11,o12,o13,o14,o15,o16;
	int ox1oy1;
	int ox1oy2,ox2oy1;
	int ox1oy3,ox3oy1;
	int ox2oy2;
	int ox2oy3,ox3oy2;
	int ox3oy3;
	double xoffmax,xoffmin,yoffmax,yoffmin,difx,dify;
	double flox,floy,fhix,fhiy,fnorm;

// intialize counters to see how pixels are overlapping
// the dimensions of the galaxy images ought to map to individual pixels
// or to an overlap of four at most so as not to waste galaxy detail

	yoffmax=-100.0;
	yoffmin=100.0;
	xoffmax=-100.0;
	xoffmin=100.0;
	oxmax=0;
	oymax=0;
	oxmin=3;
	oymin=3;
	oxoymin=10;
	oxoymax=0;
	o1=0; o2=0; o3=0; o4=0; o5=0; o6=0; o7=0; o8=0; o9=0;
	o10=0; o11=0; o12=0; o13=0; o14=0; o15=0; o16=0;
	ox1oy1=0;
	ox1oy2=0; ox2oy1=0;
	ox1oy3=0; ox3oy1=0;
	ox2oy2=0;
	ox2oy3=0; ox3oy2=0;
	ox3oy3=0;

// open the input image to be mapped onto output image

	printf("open input image: %s\n",flatimage);

	fpin=fopen(flatimage,"rb");
	if (fpin==NULL) {
		printf("Cannot open input image file: %s\n",flatimage);
		exit(1);
	}
	
	printf("input image opened correctly\n");
	printf("find corners of image mapped to picture\n");
	printf("position of UL in flat space is x:%7.3f  y:%7.3f\n",fulx,fuly);
	printf("size of image in flat space is: width:%7.3f  height%7.3f\n",
		(double)ifx*delta,(double)ify*delta);

	mul=map(fulx,fuly);
	mur=map(fulx+(double)ifx*delta,fuly);
	mll=map(fulx,fuly-(double)ify*delta);
	mlr=map(fulx+(double)ifx*delta,fuly-(double)ify*delta);

	printf("UL corner in pixels:  x=%8.5f  y=%8.5f\n",mul.ix,mul.iy);
	printf("UR corner in pixels:  x=%8.5f  y=%8.5f\n",mur.ix,mur.iy);
	printf("LL corner in pixels:  x=%8.5f  y=%8.5f\n",mll.ix,mll.iy);
	printf("LR corner in pixels:  x=%8.5f  y=%8.5f\n",mlr.ix,mlr.iy);
	printf("\n");

	printf("reading image\n");

	for(iy=0;iy<ify;iy++) {
		printf("%5d\b\b\b\b\b",iy);
		for(ix=0;ix<ifx;ix++) {
			fread(&ucval,sizeof(unsigned char),1,fpin);
			R=(double)ucval;
			fread(&ucval,sizeof(unsigned char),1,fpin);
			G=(double)ucval;
			fread(&ucval,sizeof(unsigned char),1,fpin);
			B=(double)ucval;


//-----------------------------------------------------------------------------
// map pixel

	dix=(double)ix;
	diy=(double)iy;
	ulx=fulx + dix*delta;
	uly=fuly - diy*delta;
	urx=ulx+delta;
	ury=uly;
	llx=ulx;
	lly=uly-delta;
	lrx=ulx+delta;
	lry=uly-delta;

	mul=map(ulx,uly);
	mur=map(urx,ury);
	mll=map(llx,lly);
	mlr=map(lrx,lry);

//	printf("ix:%3d  iy:%3d  R:%3.0f  G:%3.0f  B:%3.0f\n",ix,iy,R,G,B);
//	printf("upper left :  x=%8.5f  y=%8.5f\n",mul.ix,mul.iy);
//	printf("upper right:  x=%8.5f  y=%8.5f\n",mur.ix,mur.iy);
//	printf("lower left :  x=%8.5f  y=%8.5f\n",mll.ix,mll.iy);
//	printf("lower right:  x=%8.5f  y=%8.5f\n",mlr.ix,mlr.iy);
//	printf("\n");
//	getchar();




	ixm=(int)mul.ix;
	iym=(int)mul.iy;

	if ((ixm<3)||(iym<3)||(ixm>XPM1-3)||(iym>YPM1-3)) continue;

	difx=(mur.x-mul.x)/Step;
	dify=(mul.y-mll.y)/Step;
	if (xoffmax<difx) xoffmax=difx;
	if (xoffmin>difx) xoffmin=difx;
	if (yoffmax<dify) yoffmax=dify;
	if (yoffmin>dify) yoffmin=dify;

	ox=(int)(ceil(mur.ix)-floor(mul.ix));
	oy=(int)(ceil(mll.iy)-floor(mul.iy));

	if (oxmax<ox) oxmax=ox;
	if (oymax<oy) oymax=oy;
	oxoy=ox*oy;
	if (oxoymax<oxoy) oxoymax=oxoy;
	if (oxoymin>oxoy) oxoymin=oxoy;

	if(oxoy==1) o1++;
	if(oxoy==2) o2++;
	if(oxoy==3) o3++;
	if(oxoy==4) o4++;
	if(oxoy==5) o5++;
	if(oxoy==6) o6++;
	if(oxoy==7) o7++;
	if(oxoy==8) o8++;
	if(oxoy==9) o9++;
	if(oxoy==10) o10++;
	if(oxoy==11) o11++;
	if(oxoy==12) o12++;
	if(oxoy==13) o13++;
	if(oxoy==14) o14++;
	if(oxoy==15) o15++;
	if(oxoy==16) o16++;

//	if (watch==1) printf("ox=%2d  oy=%2d  oxoy=%2d\n",ox,oy,oxoy);

// single pixel

	if (oxoy==1){
		ox1oy1+=1;
		ixm=(int)mul.ix;
		iym=(int)mul.iy;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=Iscale*B*Isat/(I+Isat);
	};

// two pixels

	if ((ox==1)&&(oy==2)){
		ox1oy2+=1;
		ixm=(int)mul.ix;
		iym=(int)mul.iy;
		floy=ceil(mul.iy)-mul.iy;
		fhiy=mll.iy-floor(mll.iy);
		fnorm=floy+fhiy;

		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];	
		outR[ixm][iym]+=(floy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(floy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(floy/fnorm)*Iscale*B*Isat/(I+Isat);

		iym+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhiy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhiy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhiy/fnorm)*Iscale*B*Isat/(I+Isat);
	};

	if ((ox==2)&&(oy==1)){
		ox2oy1+=1;
		ixm=(int)mul.ix;
		iym=(int)mul.iy;
		flox=ceil(mul.ix)-mul.ix;
		fhix=mur.ix-floor(mur.ix);
		fnorm=flox+fhix;

		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(flox/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(flox/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(flox/fnorm)*Iscale*B*Isat/(I+Isat);
		ixm+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhix/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhix/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhix/fnorm)*Iscale*B*Isat/(I+Isat);
	};


// three pixels

	if ((ox==1)&&(oy==3)){
		ox1oy3+=1;
		ixm=(int)mul.ix;
		iym=(int)mul.iy;
		floy=ceil(mul.iy)-mul.iy;
		fhiy=mll.iy-floor(mll.iy);
		fnorm=1.0+floy+fhiy;

		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(floy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(floy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(floy/fnorm)*Iscale*B*Isat/(I+Isat);

		iym+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(1.0/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(1.0/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(1.0/fnorm)*Iscale*B*Isat/(I+Isat);

		iym+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhiy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhiy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhiy/fnorm)*Iscale*B*Isat/(I+Isat);
	};

	if ((ox==3)&&(oy==1)){
		ox3oy1+=1;
		ixm=(int)mul.ix;
		iym=(int)mul.iy;
		flox=ceil(mul.ix)-mul.ix;
		fhix=mur.ix-floor(mur.ix);
		fnorm=1.0+flox+fhix;

		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(flox/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(flox/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(flox/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(1.0/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(1.0/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(1.0/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhix/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhix/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhix/fnorm)*Iscale*B*Isat/(I+Isat);
	};

// four pixels

	if ((ox==2)&&(oy==2)) {
		ox2oy2+=1;
		ixm=(int)mul.ix;
		iym=(int)mul.iy;
		flox=ceil(mul.ix)-mul.ix;
		fhix=mur.ix-floor(mur.ix);
		floy=ceil(mul.iy)-mul.iy;
		fhiy=mll.iy-floor(mll.iy);
		fnorm=flox*floy + fhix*floy + flox*fhiy + fhix*fhiy;

		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(flox*floy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(flox*floy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(flox*floy/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhix*floy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhix*floy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhix*floy/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm-=1;	iym+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(flox*fhiy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(flox*fhiy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(flox*fhiy/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhix*fhiy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhix*fhiy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhix*fhiy/fnorm)*Iscale*B*Isat/(I+Isat);
	};

// six pixels

	if ((ox==2)&&(oy==3)) {
		ox2oy3+=1;
		ixm=(int)mul.ix;
		iym=(int)mul.iy;
		flox=ceil(mul.ix)-mul.ix;
		fhix=mur.ix-floor(mur.ix);
		floy=ceil(mul.iy)-mul.iy;
		fhiy=mll.iy-floor(mll.iy);
		fnorm=flox*floy+fhix*floy+flox*1.0+fhix*1.0+flox*fhiy+fhix*fhiy;

		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(flox*floy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(flox*floy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(flox*floy/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhix*floy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhix*floy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhix*floy/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm-=1; iym+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(flox/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(flox/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(flox/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhix/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhix/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhix/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm-=1; iym+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(flox*fhiy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(flox*fhiy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(flox*fhiy/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhix*fhiy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhix*fhiy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhix*fhiy/fnorm)*Iscale*B*Isat/(I+Isat);
	};	

	if ((ox==3)&&(oy==2)) {
		ox3oy2+=1;
		ixm=(int)mul.ix;
		iym=(int)mul.iy;
		flox=ceil(mul.ix)-mul.ix;
		fhix=mur.ix-floor(mur.ix);
		floy=ceil(mul.iy)-mul.iy;
		fhiy=mll.iy-floor(mll.iy);
		fnorm=flox*floy+fhix*floy+floy*1.0+fhiy*1.0+flox*fhiy+fhix*fhiy;

		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(flox*floy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(flox*floy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(flox*floy/fnorm)*Iscale*B*Isat/(I+Isat);

		iym+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(flox*fhiy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(flox*fhiy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(flox*fhiy/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1; iym-=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(floy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(floy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(floy/fnorm)*Iscale*B*Isat/(I+Isat);

		iym+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhiy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhiy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhiy/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1; iym-=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhix*floy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhix*floy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhix*floy/fnorm)*Iscale*B*Isat/(I+Isat);

		iym+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhix*fhiy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhix*fhiy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhix*fhiy/fnorm)*Iscale*B*Isat/(I+Isat);
	};	

// nine pixels

	if ((ox==3)&&(oy==3)) {

		ox3oy3+=1;
		ixm=(int)mul.ix;
		iym=(int)mul.iy;
		flox=ceil(mul.ix)-mul.ix;
		fhix=mur.ix-floor(mur.ix);
		floy=ceil(mul.iy)-mul.iy;
		fhiy=mll.iy-floor(mll.iy);
		fnorm=  flox*floy + floy + fhix*floy;
		fnorm+= flox      + 1.0  + fhix;
		fnorm+= flox*fhiy + fhiy + fhix*fhiy;

		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(flox*floy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(flox*floy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(flox*floy/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(floy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(floy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(floy/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhix*floy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhix*floy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhix*floy/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm-=2; iym+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(flox/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(flox/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(flox/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(1.0/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(1.0/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(1.0/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhix/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhix/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhix/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm-=2; iym+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(flox*fhiy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(flox*fhiy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(flox*fhiy/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhiy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhiy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhiy/fnorm)*Iscale*B*Isat/(I+Isat);

		ixm+=1;
		I=outR[ixm][iym]+outG[ixm][iym]+outB[ixm][iym];
		outR[ixm][iym]+=(fhix*fhiy/fnorm)*Iscale*R*Isat/(I+Isat);
		outG[ixm][iym]+=(fhix*fhiy/fnorm)*Iscale*G*Isat/(I+Isat);
		outB[ixm][iym]+=(fhix*fhiy/fnorm)*Iscale*B*Isat/(I+Isat);

	}; //end of nine pixels


//-----------------------------------------------------------------------------
		
		} // end of ix loop
	} // end of iy loop

	fclose(fpin);

	printf("\ngalaxy mapping statistics:\n");
	printf("xoffmax=%9.6f  xoffmin=%9.6f\n",xoffmax,xoffmin);
	printf("yoffmax=%9.6f  yoffmin=%9.6f\n",yoffmax,yoffmin);
	printf("oxmax=%4d  oymax=%4d\n",oxmax,oymax);
	printf("oxoymax=%4d  oxoymin=%4d\n",oxoymax,oxoymin);
	printf("o1=%10d\n",o1);
	printf("o2=%10d\n",o2);
	printf("o3=%10d\n",o3);
	printf("o4=%10d\n",o4);
	printf("o5=%10d\n",o5);
	printf("o6=%10d\n",o6);
	printf("o7=%10d\n",o7);
	printf("o8=%10d\n",o8);
	printf("o9=%10d\n",o9);
	printf("o10=%10d\n",o10);
	printf("o11=%10d\n",o11);
	printf("o12=%10d\n",o12);
	printf("o13=%10d\n",o13);
	printf("o14=%10d\n",o14);
	printf("o15=%10d\n",o15);
	printf("o16=%10d\n",o16);

	printf("ox1oy1=%10d\n",ox1oy1);
	printf("ox1oy2=%10d   ox2oy1=%10d sum=%10d\n",ox1oy2,ox2oy1,ox1oy2+ox2oy1);
	printf("ox2oy2=%10d\n",ox2oy2);
	printf("ox2oy3=%10d   ox3oy2=%10d sum=%10d\n",ox2oy3,ox3oy2,ox2oy3+ox3oy2);
	printf("ox3oy3=%10d\n",ox3oy3);

	printf("galaxy image mapping complete\n");


}
//*****************************************************************************
// allocation routines modified from
// "Numerical Recipes in C" by Press, Teukolsky, Vetterling and Flannery
// Cambridge University Press
// everyone who loves computing should own this book
//
/* dmatrix.c                                                           */
/*                                                                     */
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */

double **dmatrix(long nrl, long nrh, long ncl, long nch)
{
	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
	double **m;

	/* allocate pointers to rows */
	m=(double **) malloc((size_t)((nrow+1)*sizeof(double*)));
	if (!m) nrerror("allocation failure 1 in dmatrix()");
	m += 1;
	m -= nrl;

	/* allocate rows and set pointers to them */
	m[nrl]=(double *) malloc((size_t)((nrow*ncol+1)*sizeof(double)));
	if (!m[nrl]) nrerror("allocation failure 2 in dmatrix()");
	m[nrl] += 1;
	m[nrl] -= ncl;

	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

	/* return pointer to array of pointers to rows */
	return m;
}

// free a double matrix allocated by dmatrix()

void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)

{
	free((char*) (m[nrl]+ncl-1));
	free((char*) (m+nrl-1));
}
//-----------------------------------------------------------------------------
/* nrerror.c                                                         */
/*                                                                   */
/* prints an error message to the standard error handler then exits  */
/* extracted from Numerical Recipes nrutil.c                         */

void nrerror(char error_text[]) {

	char ch;

	fprintf(stderr,"Numerical Recipes run-time error...\n");
	fprintf(stderr,"%s\n",error_text);
	fprintf(stderr,"hit key to exit to system\n");
	ch=getchar();
	fprintf(stderr,"...now exiting to system...\n");
	exit(1);
}

//-----------------------------------------------------------------------------