// earth.cc
// Black Hole Earth Image Creator
// Copyright (C) 2003 Yukinori NAGATANI
// yukinori.nagatani@weizmann.ac.il

char * szTitleMessage =
	"earth: Black Hole Earth Image Creator  by Y.Nagatani";

#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <iomanip.h>
#include <fstream.h>
#include <math.h>

typedef	double	Real;
typedef	int	Bool;
#define	False	0
#define	True	1

////// PPM Raw Image Class

class PPMImage {
public:
	int		xSize,	ySize;
	unsigned char	*im;
	void allocate (int x, int y) {
		if (im != NULL)
			delete im;
		if (x < 1 || y < 1) {
			cerr << "PPMImage: initial size error as ("
			     << x << "," << y << ")." << endl;
			exit (1);
		}
		xSize = x;
		ySize = y;
		im = new unsigned char [xSize * ySize * 3];
		if (im == NULL) {
			cerr << "PPMImage: Cannot Allocate Memory."
			     << endl;
			exit (0);
		}
	}
	void clear (void) {
		unsigned char	*p = im;
		for (int n = xSize * ySize *3  ; n ; n--)
			*(p++) = 0;
	}
	PPMImage (void) {	// default constructor
		xSize = 0;
		ySize = 0;
		im = NULL;
	}
	PPMImage (int x, int y) {
		im = NULL;
		allocate (x, y);
		clear ();
	}
	PPMImage (char *szFileName) {
		xSize = 0;
		ySize = 0;
		im = NULL;
		
		ifstream	sin (szFileName);
		sin >> *this;
		sin.close ();
		
		return;
		
	}
	PPMImage (const PPMImage & ppm)	// copy constructor
	{
		im = NULL;
		allocate (ppm.xSize, ppm.ySize);
		
		unsigned char	*p = im,
				*q = ppm.im;
		for (int n = xSize * ySize *3  ; n ; n--)
			*(p++) = *(q++);
	}
	PPMImage & operator = (const PPMImage & ppm)	// substitution
	{
		if (this == &ppm)
			return	*this;
		
		im = NULL;
		allocate (ppm.xSize, ppm.ySize);
		
		unsigned char	*p = im,
				*q = ppm.im;
		for (int n = xSize * ySize *3  ; n ; n--)
			*(p++) = *(q++);
		
		return	*this;
	}
	~PPMImage (void)	// destructor
	{
		xSize = 0;
		ySize = 0;
		if (im != NULL)
			delete im;
		im = NULL;
	}
	
	Real getValue (int x, int y, int c) const {
		if (0 <= x && x < xSize &&
		    0 <= y && y < ySize) {
			y = (ySize - 1) - y;
			return	(Real) im[3*(x + xSize*y) + c] / 255;
		} else {
			return	0.0;
		}
	}
	Real getValue (Real rX, Real rY, int c) const {
		return	getValue ((int)((Real) xSize * rX),
				  (int)((Real) ySize * rY),
				  c);
	}
	PPMImage & setValue (int x, int y, int c, unsigned char val) {
		if (0 <= x && x < xSize &&
		    0 <= y && y < ySize) {
			y = (ySize - 1) - y;
			im [3*(x + xSize*y) + c] = val;
		}
		return	*this;
	}
	PPMImage & setValue (int x, int y, int c, Real val) {
		val = (val < 0.0 ? 0.0 : val);
		val = (val > 1.0 ? 1.0 : val);
		return	setValue (x, y, c, (unsigned char) (255.999 * val) );
	}
	PPMImage & setValue (Real x, Real y, int c, Real val) {
		return	setValue ((int)(((Real) xSize - 0.001) * x),
				  (int)(((Real) ySize - 0.001) * y),
				  c, val);
	}
	const PPMImage & show (void) const {
		for (int y = 0 ; y < ySize ; y++) {
			for (int x = 0 ; x < xSize ; x++) {
				cerr << "[" << (int) im[3*(x + xSize*y) + 0]
				     << "," << (int) im[3*(x + xSize*y) + 1]
				     << "," << (int) im[3*(x + xSize*y) + 2]
				     << "] ";
			}
		}
		return	*this;
	}
	friend ostream & operator << (ostream & s, const PPMImage & ppm)
	{
		s << "P6" << endl;
		s << "# Image Generated by Y.Nagatani Earth/BH" << endl;
		s << ppm.xSize << " " << ppm.ySize << endl;
		s << 255 << endl;
		
		unsigned char *p = ppm.im;
		for (int n = ppm.xSize * ppm.ySize * 3 ; n ; n--) {
			s << *(p++);
		}
		return	s;
	}
	friend istream & operator >> (istream & sin, PPMImage & ppm)
	{
		const int	nBufferSize = 256;
		char	szBuffer [nBufferSize + 1];
		int	xS,	yS;
		
		sin.getline (szBuffer, nBufferSize - 1);
		// cerr << "[" << szBuffer << "]" << endl;
		if (strcmp (szBuffer, "P6")) {
			// unmatch
			cerr << "Non PPM raw File" << endl;
			exit (1);
		}
		for (;;) {
			sin.getline (szBuffer, nBufferSize - 1);
			// cerr << "[" << szBuffer << "]" << endl;
			if (sin.eof ()) {
				cerr << "PPMImage: PPMFile Line Shortage."
				     << endl;
				exit (1);
			}
			if (szBuffer[0] == '#')
				continue;
			if (2 == sscanf (szBuffer, "%d%d", &xS, &yS)) {
				cerr << "Input Image Size is "
				     << xS << "x" << yS << endl;
				break;
			} else {
				cerr << "Cannot get Image Size."
				     << szBuffer
				     << endl;
				exit (1);
			}
		}
		sin.getline (szBuffer, nBufferSize - 1);
		if (strcmp (szBuffer, "255")) {
			// unmatch
			cerr << "Non PPM raw File with 255 depth." << endl;
			exit (1);
		}
		ppm.allocate (xS, yS);
		ppm.clear ();

		unsigned char	*p = ppm.im;
		for (int n = ppm.xSize * ppm.ySize * 3 ; n ; n--) {
			unsigned char	ucInput;
			sin.get	(ucInput);
			if (sin.eof ()) {
				cerr << "Non PPM File (Line Shortage)."
				     << endl;
				exit (1);
			}
			*(p++) = ucInput;
		}
		return	sin;
	}
};

void sample (void) {
	int		xSize = 200,	ySize = 200;
	PPMImage	imOut (xSize, ySize);
	
	for (int nX = 0 ; nX < xSize ; nX++) {
		for (int nY = 0 ; nY < ySize ; nY++) {
			Real	x = (Real) nX / xSize * 2 - 1,
				y = (Real) nY / ySize * 2 - 1;
			Real	r = sqrt (x*x + y*y) ;
			imOut.setValue (nX,nY,0, r);
			imOut.setValue (nX,nY,1, 1-r);
			imOut.setValue (nX,nY,2, x+y);
		}
	}
	
	ofstream	fout ("test01.ppm");
	fout << imOut;
	fout.close ();
	
	PPMImage	imO2 (imOut),	imO3 = imOut;
	
	fout.open ("test02.ppm");
	fout << imO2;
	fout.close ();
	
	fout.open ("test03.ppm");
	fout << imO3;
	fout.close ();
}

void rot (Real &x, Real &y, Real angle) {
	Real	X = + cos (angle) * x  + sin (angle) * y,
		Y = - sin (angle) * x  + cos (angle) * y;
	x = X;
	y = Y;
}

class EarthImage {
public:
	PPMImage	ppm;
	EarthImage & loadImage (char * szFileName = "DayEarth2.ppm") {
		ifstream	sin (szFileName);
		sin >> ppm;
		sin.close ();
	}
	EarthImage & setTile (int xSize, int ySize) {
		ppm.allocate (xSize, ySize);
		for (int x = 0 ; x < xSize ; x++)
			for (int y = 0 ; y < ySize ; y++) {
				if ( (x+y) % 2 == 0) {
					ppm.setValue(x,y,0, 1.0);
				} else {
					ppm.setValue(x,y,2, 1.0);
				}
				if (y < ySize/2)
					ppm.setValue(x,y,1, 1.0);
			}
		return	*this;
	}
	EarthImage (void) {
		setTile (12, 6);
	}
	EarthImage (int xSize, int ySize) {
		setTile (xSize, ySize);
	}
	EarthImage (char * szFileName) {
		loadImage (szFileName);
	}
	Real getValue (Real theta, Real phi, int c) const {		
		Real	X = (M_PI + phi  )/(M_PI*2);	// from 0 to 1
		Real	Y = (M_PI - theta)/(M_PI  );	// from 0 to 1
		X = X - floor(X);
		int	x = (int) (X * ((Real) ppm.xSize - 0.001));
		int	y = (int) (Y * ((Real) ppm.ySize - 0.001));
		
		return	ppm.getValue (x,y, c);
	}
};


Real phiBH (Real d) {
	Real	tanTheta = d * sqrt (27.0/4.0);
	
	if (tanTheta < 0.327161)
		return	tanTheta;
	
	return	+ 0.9779557087387333
		+ 0.5283685024076627 * tanTheta
		- 1.002411662206105  * log ( sqrt (27.0/4.0) - tanTheta );
}

void viewEarth (const EarthImage & earth, PPMImage & imOut,
		Real rLatitude, Real rLongitude, Bool isBH = False,
		Real rMag = 1.0) {
	
	int	xSize = imOut.xSize,
		ySize = imOut.ySize;
	Real	R = (isBH ? 1 : 1.0/sqrt(27.0/4.0));
	Real	xMag = sqrt(rMag),
		yMag = rMag;
	
	for (int nX = 0 ; nX < xSize ; nX++) {
		for (int nY = 0 ; nY < ySize ; nY++) {
			Real	x = (Real) nX / xSize * 2 - 1,
				y = (Real) nY / ySize * 2 - 1;
			
			if (rMag != 1.0) {			
				x = x / xMag;
				y = (1.0 - 1.0/yMag) + y/yMag;
			}
			
			Real	r = sqrt (x*x + y*y) / R;
			if (r > 1) {
				continue;
			}
			Real	theta = (isBH ? phiBH (r) : asin (r));
			Real	phi   = (r > 0 ? atan2 (y, x) : 0);
			
			Real	X = sin (theta) * cos (phi),
				Y = sin (theta) * sin (phi),
				Z = cos (theta);
			
			rot (Y,Z, M_PI/180* (+rLatitude));
			rot (Z,X, M_PI/180* (-rLongitude));
			
			Real	s = (Y < 1 ? X/sqrt(1-Y*Y) : 0);
			s = (s >= +1 ? +1 : s);
			s = (s <= -1 ? -1 : s);
			Real	thetaE = acos (Y);
			Real	phiE   = asin (s);
			
			phiE = (Z >= 0 ? phiE : M_PI - phiE);
			
			for (int c = 0 ; c < 3 ; c++) {
				Real	val = earth.getValue (thetaE,phiE ,c);
				imOut.setValue (nX,nY, c, val);
			}
		}
	}
}

int main (int nArgc, char * pszArgv []) {
	
	cerr << szTitleMessage << endl;
	
	int	nSize = 600;
	Real	rLatitude = 32,
		rLongitude = 35;
	Bool	isBH = False,
		isLoadImage = False;
	Real	rMagnify = 1.0;
	
	nArgc--;
	pszArgv++;
	for (; nArgc > 0 ; nArgc--, pszArgv++) {
		if (!strcmp (*pszArgv, "-latitude")) {
			nArgc--;
			pszArgv++;
			if (nArgc <= 0) {
				cerr << "No Argument: -latitude ##" << endl;
				exit (1);
			}
			rLatitude = atof (*pszArgv);
			continue;
		}
		if (!strcmp (*pszArgv, "-longitude")) {
			nArgc--;
			pszArgv++;
			if (nArgc <= 0) {
				cerr << "No Argument: -longitude ##" << endl;
				exit(1);
			}
			rLongitude = atof (*pszArgv);
			continue;
		}
		if (!strcmp (*pszArgv, "-size")) {
			nArgc--;
			pszArgv++;
			if (nArgc <= 0) {
				cerr << "No Argument: -size ##" << endl;
				exit(1);
			}
			nSize = atoi (*pszArgv);
			continue;
		}
		if (!strcmp (*pszArgv, "+BH")) {
			isBH = True;
			continue;
		}
		if (!strcmp (*pszArgv, "-BH")) {
			isBH = False;
			continue;
		}
		if (!strcmp (*pszArgv, "+image")) {
			isLoadImage = True;
			continue;
		}
		if (!strcmp (*pszArgv, "-magnify")) {
			nArgc--;
			pszArgv++;
			if (nArgc <= 0) {
				cerr << "No Argument: -magnify ##" << endl;
				exit(1);
			}
			rMagnify = atof (*pszArgv);
			continue;
		}
		cerr << "earth [-latitude ##] [-longitude ##] "
			"[+BH] [-BH] [+image] [-magnify ##]" << endl;
		cerr << "Invalid Argument: " << *pszArgv << endl;
		exit (1);
	}
	
	EarthImage	earth;
	if (isLoadImage)
		earth.loadImage ("DayEarth2.ppm");
	
	PPMImage	imOut (nSize,nSize);
	
	//viewEarth (earth, imOut, 32, 35, True);
	viewEarth (earth, imOut, rLatitude, rLongitude, isBH, rMagnify);
	
	// ofstream	fout ("test.ppm");
	cout << imOut;
	// fout.close ();
	
	return 0;
}