Sign in to follow this  
SiS-Shadowman

Calculate Sun's Position from Earth POV

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 this post


Link to post
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
Admiral

Share this post


Link to post
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.

please help me :)

Share this post


Link to post
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[i][TERM_A] * cos(terms[i][TERM_B] + terms[i][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[i]*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[i] = EarthPeriodicTermSum(L_TERMS[i], l_subcount[i], 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[i] = EarthPeriodicTermSum(B_TERMS[i], b_subcount[i], 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[i] = EarthPeriodicTermSum(R_TERMS[i], r_subcount[i], 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[i][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[i][TERM_PSI_A] + jce*PE_TERMS[i][TERM_PSI_B]) * sin(xy_term_sum);
sum_epsilon += (PE_TERMS[i][TERM_EPS_C] + jce*PE_TERMS[i][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) );
}



The header file:

#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 40

enum { 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 this post


Link to post
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]

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this