{
	Coordinate mapping for point above Earth
	assuming uniform projection of latitude and longitude

		5/1/96		Greg Ward

	A1 = latitude position (degrees)
	A2 = longitude position (degrees)
	A3 = height above Earth (kilometers)
	A4 = Julian date (1-365)
	A5 = Greenwich Mean Time (0-24)

	By modifying this file so that the above arguments are defined
	instead as constants, significant gains in speed are possible.
	(Be sure to change Su, Sv and Sw to constant definitions also.)

	Y-axis points North, Z-axis points up (away from Earth)
}
{
	Example usage:

	void colorpict epict
	7 Eval Eval Eval Earth.pic Earth.cal u v
	0
	5 lat long height JDate GMT

	epict glow eglow
	0
	0
	4 1 1 1 0

	eglow source earth
	0
	0
	4 0 0 -1 max_angle

	where:
		max_angle = 2*180/PI * asin(EarthRad/(EarthRad+height))

	If you are really close to the Earth (below atmosphere), use:

	!gensky month day time -a lat -o lon -m std_meridian

	skybright colorpict epict
	7 noop noop noop Earth.pic Earth.cal u v
	0
	3 lat long height

	epict glow eglow
	0
	0
	4 1 1 1 0

	eglow source earth
	0
	0
	4 0 0 -1 max_angle
}

EarthRad : 6378.5;				{ Earth radius in kilometers }
EarthTilt : 23.45*DEGREE;			{ Earth's orbital axis tilt }
EarthIrrad : 770 {watts/meter^2};		{ Maximum Earth irradiance }
PA : 256/512;					{ Picture Height/Width }

		{ Change the following to constants to speed things up }
Lat = arg(1)*DEGREE;		{ Position latitude }
Lon = arg(2)*DEGREE;		{ Position longitude }
Height = arg(3);		{ Position height above Earth (km) }
JDate = arg(4);			{ Julian date (1-365) }
GMT = arg(5);			{ Greenwich Mean Time (0-24) }

GOu = EarthRad*sin(Lon)*cos(Lat);		{ Starting Earth position }
GOv = EarthRad*cos(Lon)*cos(Lat);
GOw = EarthRad*sin(Lat);

raduv = max(FTINY, sqrt(GOu*GOu + GOv*GOv));	{ Earthwise direction vector }
GVXu = -GOv / raduv;
GVXv = GOu / raduv;
GVXw = 0;
GVYu = (-GOw*GOu) / EarthRad / raduv;
GVYv = (GOw*-GOv) / EarthRad / raduv;
GVYw = (GOu*GOu + GOv*GOv) / EarthRad / raduv;
radxy = max(FTINY, sqrt(Dx*Dx + Dy*Dy));
GDu = (GVXu*Dx + GVYu*Dy) / radxy;
GDv = (GVXv*Dx + GVYv*Dy) / radxy;
GDw = (GVXw*Dx + GVYw*Dy) / radxy;

tana = if(FTINY+Dz, radxy/FTINY, radxy/-Dz);

det = sq(EarthRad+Height) - (tana*tana + 1)*(2*EarthRad*Height + Height*Height);

a = if(tana-FTINY, b/tana, Height);

b = if(tana-FTINY, (EarthRad + Height - Sqrt(det)) / (tana + 1/tana), 0);

c = EarthRad + Height - a;

digfactor = 1 - (a - Height)/EarthRad;
						{ Global position vector }
GPu = GOu*digfactor + GDu*b;
GPv = GOv*digfactor + GDv*b;
GPw = GOw*digfactor + GDw*b;
						{ Earthly coordinates }
longitude = atan2(GPu, GPv);
latitude = Asin(GPw/EarthRad);
						{ Image lookup coordinates }
u = if(PA-1, 1, 1/PA) * (.5 - longitude/(2*PI));
v = if(PA-1, PA, 1) * (.5 + latitude/PI);
						{ Solar declination (lat.) }
Slat = EarthTilt * sin(2*PI/368 * (JDate - 81));
						{ Solar longitude }
Slon = 2*PI/24 * (GMT - 12);
						{ Solar position vector }
Su = EarthRad*sin(Slon)*cos(Slat);
Sv = EarthRad*cos(Slon)*cos(Slat);
Sw = EarthRad*sin(Slat);
						{ Solar cosine }
Scos = (Su*GPu + Sv*GPv + Sw*GPw)/(EarthRad*EarthRad);
						{ Earth radiance multiplier }
Eradmult = EarthIrrad / PI * noneg(Scos);
						{ Earth radiance }
Eval(rho) = Eradmult * rho;