# Calculate Sun's Position from Earth POV

This topic is 4454 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

okay, i'm back with my little game and back with my little problems :( i need to calculate the sun's position for an earth viewer, but my math isn't good enough for that -.- the orbit from earth should just look like a circle, but it is turned by the latitude of the viewer. how can i accomplish that? can somebody please post a formula to do this?

##### Share on other sites
How accurate do you need this to be? It should be very easy to create an oblique (~23°) vector around which to uniformly rotate the sun using your APIs matrix functions. If you need to model the sun for seasons and so forth then there are a few more considerations to accommodate.

Regards

##### Share on other sites
this would be more than ok. my current function is really incorrect ( the sun doesn't move in a circle but in some parabola shape over the sky *fg )

i don't need any model. i just need a formula to calculate the sun's position and perhaps its rotation. the sun itself is just a colored quad, the sunshine is rendered on the skydome. to do this, i need a vector with its position.

##### Share on other sites
Here's something with a "little" more precision :) Might be a little overkill though ;)

It uses an external library for vector math, but it shoud be straightforward to replace the few method calls by your own.

static double LimitDegrees(double degrees){	degrees /= 360.0;	double limited = 360.0 * (degrees-floor(degrees));	if( limited < 0 ) limited += 360.0;		return( limited );}static double ThirdOrderPolynomial(double a, double b, double c, double d, double x){	return( ((a*x + b)*x + c)*x + d );}static double JulianDay(int year, int month, int day, int hour, int minute, int second, double tz){	double day_decimal = day + (hour - tz + (minute + second / 60.0) / 60.0) / 24.0;		if( month < 3 ) {		month += 12;		year--;	}		double julian_day = floor(365.25 * (year + 4716.0)) + floor(30.6001 * (month + 1)) + day_decimal - 1524.5;		if( julian_day > 2299160.0 ) {		int a = (int)floor((float)(year / 100));		julian_day += (2 - a + (int)floor((float)(a/4)));	}	 	return( julian_day );}static double EarthPeriodicTermSum(const double terms[][TERM_COUNT], int count, double jme){	double sum = 0;		for( int i = 0; i < count; i++ ) sum += terms[TERM_A] * cos(terms[TERM_B] + terms[TERM_C] * jme);		return( sum );}		static double EarthValues(double term_sum[], int count, double jme){	double sum = 0;		for( int i = 0; i < count; i++ ) sum += term_sum*pow(jme, i);	return( sum / 1.0e8 );}static double EarthHeliocentricLongitude(double jme){	double sum[L_COUNT];		for( int i = 0; i < L_COUNT; i++ ) sum = EarthPeriodicTermSum(L_TERMS, l_subcount, jme);			return( LimitDegrees( EarthValues(sum, L_COUNT, jme) * sgcn<double>::rad2deg() ) );}static double EarthHeliocentricLatitude(double jme){	double sum[B_COUNT];		for( int i = 0; i < B_COUNT; i++ ) sum = EarthPeriodicTermSum(B_TERMS, b_subcount, jme);			return( EarthValues(sum, B_COUNT, jme) * sgcn<double>::rad2deg() );}static double EarthRadiusVector(double jme){	double sum[R_COUNT];		for( int i = 0; i < R_COUNT; i++ ) sum = EarthPeriodicTermSum(R_TERMS, r_subcount, jme);			return( EarthValues(sum, R_COUNT, jme) );}static double XYTermSum(int i, double x[TERM_X_COUNT]){	double sum = 0;		for( int j = 0; j < TERM_Y_COUNT; j++ ) sum += x[j]*Y_TERMS[j];			return( sum );}static void NutationLongitudeObliquity(double jce, double x[TERM_X_COUNT], double *del_psi, double *del_epsilon){	double xy_term_sum, sum_psi = 0, sum_epsilon = 0;		for( int i = 0; i < Y_COUNT; i++ ) {		xy_term_sum = XYTermSum(i, x) * sgcn<double>::deg2rad();		sum_psi += (PE_TERMS[TERM_PSI_A] + jce*PE_TERMS[TERM_PSI_B]) * sin(xy_term_sum);		sum_epsilon += (PE_TERMS[TERM_EPS_C] + jce*PE_TERMS[TERM_EPS_D]) * cos(xy_term_sum);	}		*del_psi = sum_psi / 36000000.0;	*del_epsilon = sum_epsilon / 36000000.0;}static double EclipticMeanObliquity(double jme){	double u = jme / 10.0;		return( 84381.448 + u*(-4680.96 + u*(-1.55 + u*(1999.25 + u*(-51.38 + u*(-249.67 + u*(  -39.05 + u*( 7.12 + u*(  27.87 + u*(  5.79 + u*2.45))))))))) );}static double GeocentricSunRightAscension(double lamda, double epsilon, double beta){	double lamda_rad = lamda * sgcn<double>::deg2rad();	double epsilon_rad = epsilon * sgcn<double>::deg2rad();	double alpha = atan2(sin(lamda_rad)*cos(epsilon_rad) - tan(beta * sgcn<double>::deg2rad())*sin(epsilon_rad), cos(lamda_rad));		return( LimitDegrees(alpha * sgcn<double>::rad2deg()) );}  static double GeocentricSunDeclination(double beta, double epsilon, double lamda){	double beta_rad = beta * sgcn<double>::deg2rad();	double epsilon_rad = epsilon * sgcn<double>::deg2rad();	double delta = asin(sin(beta_rad)*cos(epsilon_rad) + cos(beta_rad)*sin(epsilon_rad)*sin(lamda * sgcn<double>::deg2rad()));		return( delta * sgcn<double>::rad2deg() );}static void RAParallaxAndTopocentricDec(double latitude, double elevation, double xi, double h, double delta, double *delta_alpha, double *delta_prime){	double lat_rad = latitude * sgcn<double>::deg2rad();	double xi_rad = xi * sgcn<double>::deg2rad();	double h_rad = h * sgcn<double>::deg2rad();	double delta_rad = delta * sgcn<double>::deg2rad();	double u = atan(0.99664719 * tan(lat_rad));	double y = 0.99664717 * sin(u) + elevation*sin(lat_rad)/6378140.0;	double x = cos(u) + elevation*cos(lat_rad)/6378140.0;		*delta_alpha = atan2(-x*sin(xi_rad)*sin(h_rad), cos(delta_rad) - x*sin(xi_rad)*cos(h_rad)) * sgcn<double>::rad2deg();	*delta_prime = atan2((sin(delta_rad) - y*sin(xi_rad))*cos(*delta_alpha), cos(delta_rad) - x*sin(xi_rad) *cos(h_rad)) * sgcn<double>::rad2deg();}static double TopocentricElevationAngle(double latitude, double delta_prime, double h_prime){	double lat_rad = latitude * sgcn<double>::deg2rad();	double delta_prime_rad = delta_prime * sgcn<double>::deg2rad();	double e0 = asin(sin(lat_rad)*sin(delta_prime_rad) + cos(lat_rad)*cos(delta_prime_rad) * cos(h_prime * sgcn<double>::deg2rad()));		return( e0 * sgcn<double>::rad2deg() );}static double TopocentricAzimuthAngle(double h_prime, double latitude, double delta_prime){	double h_prime_rad = h_prime * sgcn<double>::deg2rad();	double lat_rad = latitude * sgcn<double>::deg2rad();	double az = atan2(sin(h_prime_rad), cos(h_prime_rad)*sin(lat_rad) - tan(delta_prime * sgcn<double>::deg2rad())*cos(lat_rad));		return( az * sgcn<double>::rad2deg() );}fnpolar ComputeSunAngle(const SObserver *O){	double x[TERM_X_COUNT];	// Calculate Julian day and century	const DateTime &dt = O->datetime;	double jd = JulianDay(dt.year(), dt.month(), dt.day(), dt.hour(), dt.minute(), dt.second(), O->timezone);	double jc = (jd - 2451545.0) / 36525.0;	// Calculate Julian Ephemeris day, century, and millenium	double jde = jd + O->delta_t / 86400.0;	double jce = (jde - 2451545.0) / 36525.0;	double jme = jce / 10.0;	// Compute earth heliocentric longitude, latitude, and radius vector	double l = EarthHeliocentricLongitude(jme);	double b = EarthHeliocentricLatitude(jme);	double r = EarthRadiusVector(jme);	// Compute geocentric longitude and latitude	double theta = l + 180.0; if( theta >= 360.0 ) theta -= 360.0;	double beta = -b;	// Calculate the nutation in longitude and obliquity	x[TERM_X0] = ThirdOrderPolynomial(1.0 / 189474.0, -0.0019142, 445267.11148, 297.85036, jce);		// mean elongation of the moon from the sun	x[TERM_X1] = ThirdOrderPolynomial(-1.0 / 300000.0, -0.0001603, 35999.05034, 357.52772, jce);		// mean anomaly of the sun	x[TERM_X2] = ThirdOrderPolynomial(1.0 / 56250.0, 0.0086972, 477198.867398, 134.96298, jce);		// mean anomaly of the moon	x[TERM_X3] = ThirdOrderPolynomial(1.0 / 327270.0, -0.0036825, 483202.017538, 93.27191, jce);		// moon's argument of latitude	x[TERM_X4] = ThirdOrderPolynomial(1.0 / 450000.0, 0.0020708, -1934.136261, 125.04452, jce);		// longitude of the ascending node of the moon's mean orbit on the ecliptic	double delta_psi, delta_epsilon;	NutationLongitudeObliquity(jce, x, &delta_psi, &delta_epsilon);	// Calculate the true obliquity of the ecliptic	double epsilon0 = EclipticMeanObliquity(jme);	double epsilon = delta_epsilon + epsilon0 / 3600.0;	// Calculate the apparent sidereal time at Greenwich at any given time	double delta_tau = -20.4898 / (3600.0*r);	double lamda = theta + delta_psi + delta_tau;	double nu0 = LimitDegrees(280.46061837 + 360.98564736629 * (jd - 2451545.0) + jc*jc*(0.000387933 - jc/38710000.0));	double nu = nu0 + delta_psi * cos(epsilon * sgcn<double>::deg2rad());		// Calculate the geocentric sun right ascension and declination	double alpha = GeocentricSunRightAscension(lamda, epsilon, beta);	double delta = GeocentricSunDeclination(beta, epsilon, lamda);	// Calculate the observer local hour angle	double h = nu + O->longitude - alpha;	// Calculate the topocentric sun right ascension	double xi = 8.794 / (3600.0 * r);	double delta_alpha, delta_prime;	RAParallaxAndTopocentricDec(O->latitude, O->altitude, xi, h, delta, &delta_alpha, &delta_prime);	double alpha_prime = alpha + delta_alpha;	// Calculate the topocentric zenith angle	double h_prime = h - delta_alpha;	double e = TopocentricElevationAngle(O->latitude, delta_prime, h_prime);	// TODO: apply refractive correction at standard air pressure and temperature here ?	double zenith = 90.0 - e;	// Calculate the topocentric azimuth angle	double azimuth = TopocentricAzimuthAngle(h_prime, O->latitude, delta_prime) + 180.0;	// Return angles as normalized polar coordinates	return( fnpolar((float)azimuth, (float)zenith) );}

#define L_COUNT		6#define B_COUNT		2#define R_COUNT		5#define Y_COUNT		63#define L_MAX_SUBCOUNT		64#define B_MAX_SUBCOUNT		5#define R_MAX_SUBCOUNT		40enum { TERM_A = 0, TERM_B = 1, TERM_C = 2, TERM_COUNT = 3 };enum { TERM_X0 = 0, TERM_X1 = 1, TERM_X2 = 2, TERM_X3 = 3, TERM_X4 = 4, TERM_X_COUNT = 5, TERM_Y_COUNT = 5 };enum { TERM_PSI_A = 0, TERM_PSI_B = 1, TERM_EPS_C = 2, TERM_EPS_D = 3, TERM_PE_COUNT = 4 };static const double L_TERMS[L_COUNT][L_MAX_SUBCOUNT][TERM_COUNT] ={	{		{175347046.0,0,0},		{3341656.0,4.6692568,6283.07585},		{34894.0,4.6261,12566.1517},		{3497.0,2.7441,5753.3849},		{3418.0,2.8289,3.5231},		{3136.0,3.6277,77713.7715},		{2676.0,4.4181,7860.4194},		{2343.0,6.1352,3930.2097},		{1324.0,0.7425,11506.7698},		{1273.0,2.0371,529.691},		{1199.0,1.1096,1577.3435},		{990,5.233,5884.927},		{902,2.045,26.298},		{857,3.508,398.149},		{780,1.179,5223.694},		{753,2.533,5507.553},		{505,4.583,18849.228},		{492,4.205,775.523},		{357,2.92,0.067},		{317,5.849,11790.629},		{284,1.899,796.298},		{271,0.315,10977.079},		{243,0.345,5486.778},		{206,4.806,2544.314},		{205,1.869,5573.143},		{202,2.4458,6069.777},		{156,0.833,213.299},		{132,3.411,2942.463},		{126,1.083,20.775},		{115,0.645,0.98},		{103,0.636,4694.003},		{102,0.976,15720.839},		{102,4.267,7.114},		{99,6.21,2146.17},		{98,0.68,155.42},		{86,5.98,161000.69},		{85,1.3,6275.96},		{85,3.67,71430.7},		{80,1.81,17260.15},		{79,3.04,12036.46},		{71,1.76,5088.63},		{74,3.5,3154.69},		{74,4.68,801.82},		{70,0.83,9437.76},		{62,3.98,8827.39},		{61,1.82,7084.9},		{57,2.78,6286.6},		{56,4.39,14143.5},		{56,3.47,6279.55},		{52,0.19,12139.55},		{52,1.33,1748.02},		{51,0.28,5856.48},		{49,0.49,1194.45},		{41,5.37,8429.24},		{41,2.4,19651.05},		{39,6.17,10447.39},		{37,6.04,10213.29},		{37,2.57,1059.38},		{36,1.71,2352.87},		{36,1.78,6812.77},		{33,0.59,17789.85},		{30,0.44,83996.85},		{30,2.74,1349.87},		{25,3.16,4690.48}	},	{		{628331966747.0,0,0},		{206059.0,2.678235,6283.07585},		{4303.0,2.6351,12566.1517},		{425.0,1.59,3.523},		{119.0,5.796,26.298},		{109.0,2.966,1577.344},		{93,2.59,18849.23},		{72,1.14,529.69},		{68,1.87,398.15},		{67,4.41,5507.55},		{59,2.89,5223.69},		{56,2.17,155.42},		{45,0.4,796.3},		{36,0.47,775.52},		{29,2.65,7.11},		{21,5.34,0.98},		{19,1.85,5486.78},		{19,4.97,213.3},		{17,2.99,6275.96},		{16,0.03,2544.31},		{16,1.43,2146.17},		{15,1.21,10977.08},		{12,2.83,1748.02},		{12,3.26,5088.63},		{12,5.27,1194.45},		{12,2.08,4694},		{11,0.77,553.57},		{10,1.3,3286.6},		{10,4.24,1349.87},		{9,2.7,242.73},		{9,5.64,951.72},		{8,5.3,2352.87},		{6,2.65,9437.76},		{6,4.67,4690.48}	},	{		{52919.0,0,0},		{8720.0,1.0721,6283.0758},		{309.0,0.867,12566.152},		{27,0.05,3.52},		{16,5.19,26.3},		{16,3.68,155.42},		{10,0.76,18849.23},		{9,2.06,77713.77},		{7,0.83,775.52},		{5,4.66,1577.34},		{4,1.03,7.11},		{4,3.44,5573.14},		{3,5.14,796.3},		{3,6.05,5507.55},		{3,1.19,242.73},		{3,6.12,529.69},		{3,0.31,398.15},		{3,2.28,553.57},		{2,4.38,5223.69},		{2,3.75,0.98}	},	{		{289.0,5.844,6283.076},		{35,0,0},		{17,5.49,12566.15},		{3,5.2,155.42},		{1,4.72,3.52},		{1,5.3,18849.23},		{1,5.97,242.73}	},	{		{114.0,3.142,0},		{8,4.13,6283.08},		{1,3.84,12566.15}	},	{		{1,3.14,0}	}};static const double B_TERMS[B_COUNT][B_MAX_SUBCOUNT][TERM_COUNT] ={	{		{280.0,3.199,84334.662},		{102.0,5.422,5507.553},		{80,3.88,5223.69},		{44,3.7,2352.87},		{32,4,1577.34}	},	{		{9,3.9,5507.55},		{6,1.73,5223.69}	}};static const double R_TERMS[R_COUNT][R_MAX_SUBCOUNT][TERM_COUNT] ={	{		{100013989.0,0,0},		{1670700.0,3.0984635,6283.07585},		{13956.0,3.05525,12566.1517},		{3084.0,5.1985,77713.7715},		{1628.0,1.1739,5753.3849},		{1576.0,2.8469,7860.4194},		{925.0,5.453,11506.77},		{542.0,4.564,3930.21},		{472.0,3.661,5884.927},		{346.0,0.964,5507.553},		{329.0,5.9,5223.694},		{307.0,0.299,5573.143},		{243.0,4.273,11790.629},		{212.0,5.847,1577.344},		{186.0,5.022,10977.079},		{175.0,3.012,18849.228},		{110.0,5.055,5486.778},		{98,0.89,6069.78},		{86,5.69,15720.84},		{86,1.27,161000.69},		{85,0.27,17260.15},		{63,0.92,529.69},		{57,2.01,83996.85},		{56,5.24,71430.7},		{49,3.25,2544.31},		{47,2.58,775.52},		{45,5.54,9437.76},		{43,6.01,6275.96},		{39,5.36,4694},		{38,2.39,8827.39},		{37,0.83,19651.05},		{37,4.9,12139.55},		{36,1.67,12036.46},		{35,1.84,2942.46},		{33,0.24,7084.9},		{32,0.18,5088.63},		{32,1.78,398.15},		{28,1.21,6286.6},		{28,1.9,6279.55},		{26,4.59,10447.39}	},	{		{103019.0,1.10749,6283.07585},		{1721.0,1.0644,12566.1517},		{702.0,3.142,0},		{32,1.02,18849.23},		{31,2.84,5507.55},		{25,1.32,5223.69},		{18,1.42,1577.34},		{10,5.91,10977.08},		{9,1.42,6275.96},		{9,0.27,5486.78}	},	{		{4359.0,5.7846,6283.0758},		{124.0,5.579,12566.152},		{12,3.14,0},		{9,3.63,77713.77},		{6,1.87,5573.14},		{3,5.47,18849}	},	{		{145.0,4.273,6283.076},		{7,3.92,12566.15}	},	{		{4,2.56,6283.08}	}};static const int Y_TERMS[Y_COUNT][TERM_Y_COUNT] ={	{0,0,0,0,1},	{-2,0,0,2,2},	{0,0,0,2,2},	{0,0,0,0,2},	{0,1,0,0,0},	{0,0,1,0,0},	{-2,1,0,2,2},	{0,0,0,2,1},	{0,0,1,2,2},	{-2,-1,0,2,2},	{-2,0,1,0,0},	{-2,0,0,2,1},	{0,0,-1,2,2},	{2,0,0,0,0},	{0,0,1,0,1},	{2,0,-1,2,2},	{0,0,-1,0,1},	{0,0,1,2,1},	{-2,0,2,0,0},	{0,0,-2,2,1},	{2,0,0,2,2},	{0,0,2,2,2},	{0,0,2,0,0},	{-2,0,1,2,2},	{0,0,0,2,0},	{-2,0,0,2,0},	{0,0,-1,2,1},	{0,2,0,0,0},	{2,0,-1,0,1},	{-2,2,0,2,2},	{0,1,0,0,1},	{-2,0,1,0,1},	{0,-1,0,0,1},	{0,0,2,-2,0},	{2,0,-1,2,1},	{2,0,1,2,2},	{0,1,0,2,2},	{-2,1,1,0,0},	{0,-1,0,2,2},	{2,0,0,2,1},	{2,0,1,0,0},	{-2,0,2,2,2},	{-2,0,1,2,1},	{2,0,-2,0,1},	{2,0,0,0,1},	{0,-1,1,0,0},	{-2,-1,0,2,1},	{-2,0,0,0,1},	{0,0,2,2,1},	{-2,0,2,0,1},	{-2,1,0,2,1},	{0,0,1,-2,0},	{-1,0,1,0,0},	{-2,1,0,0,0},	{1,0,0,0,0},	{0,0,1,2,0},	{0,0,-2,2,2},	{-1,-1,1,0,0},	{0,1,1,0,0},	{0,-1,1,2,2},	{2,-1,-1,2,2},	{0,0,3,2,2},	{2,-1,0,2,2},};static const double PE_TERMS[Y_COUNT][TERM_PE_COUNT] ={	{-171996,-174.2,92025,8.9},	{-13187,-1.6,5736,-3.1},	{-2274,-0.2,977,-0.5},	{2062,0.2,-895,0.5},	{1426,-3.4,54,-0.1},	{712,0.1,-7,0},	{-517,1.2,224,-0.6},	{-386,-0.4,200,0},	{-301,0,129,-0.1},	{217,-0.5,-95,0.3},	{-158,0,0,0},	{129,0.1,-70,0},	{123,0,-53,0},	{63,0,0,0},	{63,0.1,-33,0},	{-59,0,26,0},	{-58,-0.1,32,0},	{-51,0,27,0},	{48,0,0,0},	{46,0,-24,0},	{-38,0,16,0},	{-31,0,13,0},	{29,0,0,0},	{29,0,-12,0},	{26,0,0,0},	{-22,0,0,0},	{21,0,-10,0},	{17,-0.1,0,0},	{16,0,-8,0},	{-16,0.1,7,0},	{-15,0,9,0},	{-13,0,7,0},	{-12,0,6,0},	{11,0,0,0},	{-10,0,5,0},	{-8,0,3,0},	{7,0,-3,0},	{-7,0,0,0},	{-7,0,3,0},	{-7,0,3,0},	{6,0,0,0},	{6,0,-3,0},	{6,0,-3,0},	{-6,0,3,0},	{-6,0,3,0},	{5,0,0,0},	{-5,0,3,0},	{-5,0,3,0},	{-5,0,3,0},	{4,0,0,0},	{4,0,0,0},	{4,0,0,0},	{-4,0,0,0},	{-4,0,0,0},	{-4,0,0,0},	{3,0,0,0},	{-3,0,0,0},	{-3,0,0,0},	{-3,0,0,0},	{-3,0,0,0},	{-3,0,0,0},	{-3,0,0,0},	{-3,0,0,0},};static const int l_subcount[L_COUNT]	= { 64, 34, 20, 7, 3, 1 };static const int b_subcount[B_COUNT]	= { 5, 2 };static const int r_subcount[R_COUNT]	= { 40, 10, 6, 2, 1 };

The observer metrics structure:
	struct SObserver {			DateTime		datetime;		// Observer time and date			double			delta_t;			// difference between earth rotation time and terrestrial time (from observation) (in seconds)			double			timezone;		// observer timezone (negative west of greenwich) (in hours)			double			longitude;		// observer longitude (negative west of greenwich) (in degrees)			double			latitude;		// observer latitude (negative south of equator) (in degrees)			double 			altitude;		// observer altitude (in meters)	};

##### Share on other sites
wow :O

i think it will take some time till i understand what you are doing there *g
thanks for the source :)

EDIT:

sorry, but i think that goes just beyond my understanding. i don't need anything as precise as your solution. perhaps someone can post an easier formula to calculate the position? please :)

[Edited by - SiS-Shadowman on October 8, 2006 6:59:41 AM]

1. 1
2. 2
Rutin
20
3. 3
khawk
16
4. 4
A4L
14
5. 5

• 11
• 16
• 26
• 10
• 11
• ### Forum Statistics

• Total Topics
633756
• Total Posts
3013707
×