# How to calculate the Sun's declination for a specific location based on the axial tilt of Earth throughout the year?

How to calculate the Sun's declination for a specific location (i.e. relative to a specific coordinate) based on the axial tilt of Earth throughout the year? Any algorithm for code available in this regard?

So I tried to port the below JS code to cpp

 #include <iostream>
#include <cmath>

using namespace std;

double earthRotationAngle(double jd);
double greenwichMeanSiderealTime(double jd);
double Julian_day(string date);

double pi = atan(1)*4, lat, lon;
string date;

double Julian_day(string date){

int Y = stoi(date.substr(0, 4));
int M = stoi(date.substr(5, 2));
int D = stoi(date.substr(8));

double jd = (1461 * (Y + 4800 + (M - 14)/12))/4 + (367 * (M - 2 - 12 * ((M - 14)/12)))/ 12 - (3 * ((Y + 4900 + (M - 14)/12)/100))/4 + D - 32075;
/*double A = Y/100;
double B = A/4;
double C = 2-A+B;
double E = 365.25*(Y+4716);
double F = 30.6001*(M+1);
double jd= C+D+E+F-1524.5;*/

return jd;
}

class sunPosition {
// Takes Julian date and returns Sun's right ascension and declination
public:
double n=Julian_day(date)-2451545.0;

double L=fmod((280.460+0.9856474*n),360);

double Lg(){
if(L < 0){L+=360;}
if(g < 0){g+=pi*2.0;}
return lamba;
}

//double beta = 0;
double ra = atan2(cos(eps)*sin(Lg()),cos(Lg()));
double dec = asin(sin(eps)*sin(Lg()));

double ara(){
if(ra < 0){ra+=pi*2;}

return a_ra;
}

};

//Greg Miller ([email protected]) 2021
//Released as public domain
//http://www.celestialprogramming.com/

//All input and output angles are in radians, jd is Julian Date in UTC
//The inputs are right ascension, sun's declination, latitude, longitude and Julian date
//(double ra,double dec,double lat,double lon,double jd_ut)

//Meeus 13.5 and 13.6, modified so West longitudes are negative and 0 is North
public:
const double gmst=greenwichMeanSiderealTime(Julian_day(date));
double localSiderealTime = fmod((gmst+lon),(2*pi));

double H = (localSiderealTime - ara());

double Hfix(){
if(H<0){H+=2*pi;}
if(H>pi){H=H-2*pi;}
return H;
}

double az = atan2(sin(Hfix()), (cos(Hfix())*sin(lat) - tan(a_dec)*cos(lat)));
double a = asin(sin(lat)*sin(a_dec) + cos(lat)*cos(a_dec)*cos(Hfix()));

double azfix(){
az -= pi;
if(az<0){az+=2*pi;}
return az;
}
//returns (az,a,localSiderealTime,H);
};

double greenwichMeanSiderealTime(double jd){
//"Expressions for IAU 2000 precession quantities" N. Capitaine1,P.T.Wallace2, and J. Chapront
const double t = ((jd - 2451545.0)) / 36525.0;

double gmst = earthRotationAngle(jd)+(0.014506 + 4612.156534*t + 1.3915817*pow(t,2) - 0.00000044 *pow(t,3)- 0.000029956*pow(t,4) - 0.0000000368*pow(t,5))/60.0/60.0*pi/180.0;  //eq 42
gmst=fmod(gmst,(2*pi));
if(gmst<0) gmst+=2*pi;

return gmst;
}

double earthRotationAngle(double jd){
//IERS Technical Note No. 32

const double t = jd- 2451545.0;
const double f = fmod(jd,1.0);

double theta = 2*pi * (f + 0.7790572732640 + 0.00273781191135448 * t); //eq 14
theta=fmod(theta,(2*pi));
if(theta<0){theta+=2*pi;}

return theta;
}

int main(){
cout << "Enter your location latitude:  " <<endl;
cin >> lat;
cout << "Enter your location longitude: " <<endl;
cin >> lon;
cout << "Enter the day of the year as yyyy-mm-dd: "<< endl;
cin >> date; //= "2019-01-09";

//sunPosition sun;

cout << endl <<"Julian_day: " << Julian_day(date) << endl;
cout << "earthRotationAngle: " << earthRotationAngle(Julian_day(date)) << endl;
cout << "Azimuth: " <<alt.azfix() << endl;                      // azimuth
cout << "Altitude: " <<alt.a << endl;                           // altitude may be
cout << "Local Sidereal time: " <<alt.localSiderealTime << endl;    // localSiderealTime
cout << "Hour angle: " <<alt.Hfix() << endl;                    // hour angle
cout << "Right ascension: " << alt.ara() << endl;               // Right ascension
cout << "Sun's declination: " << alt.a_dec << endl;             // Sun's declination
}



and it returned for lat 77.2069702520977 lon 118.639627806683 date: 1472-Aug-18

Julian_day: 2.25893e+06

earthRotationAngle: 2.68419

Azimuth: 1.57751

Altitude: 0.604463

Local Sidereal time: 1.82557

Hour angle: -1.72583

Right ascension: 9.83458

Sun's declination: 13.1531

Which do not coincide with this test data. So any more insight on it?

• As Greg says below, the Sun's declination only depends on the time. It doesn't depend on the observer's location. Feb 15, 2022 at 8:15
• Greg's equations are for the era 1950 to 2050. They lose accuracy outside that era. BTW, your Julian day number doesn't have enough digits. The JDN of Noon 1472-Aug-18 is 2258927 if you're using the Gregorian calandar, or 2258936 if you're using the Julian calendar. I should mention that Horizons uses the Julian calendar for dates before 1582-Oct-15. Feb 22, 2022 at 10:36

You've asked for "declination", but the rest of the question sounds like you're interested more in alt/az coordinates (declination doesn't change based on location). Since you need to compute the declination for alt/az coordinates, here is how to do both.

The date input to these functions are Julian Dates. The code is in JavaScript.

First compute the sun's right ascension and declination:

//Low precision sun position from Astronomical Almanac page C5 (2017 ed).
//Accuracy 1deg from 1950-2050
function sunPosition(jd)    {
n=jd-2451545.0;
L=(280.460+0.9856474*n)%360;
if(L<0){L+=360;}
if(g<0){g+=Math.PI*2.0;}

beta=0.0;
ra=Math.atan2(Math.cos(eps)*Math.sin(lamba),Math.cos(lamba));
dec=Math.asin(Math.sin(eps)*Math.sin(lamba));
if(ra<0){ra+=Math.PI*2;}
}


Then compute the Alt/Az position of those RA Dec coordinates for a given latitude and longitude:

//Greg Miller ([email protected]) 2021
//Released as public domain
//http://www.celestialprogramming.com/

//All input and output angles are in radians, jd is Julian Date in UTC
//Meeus 13.5 and 13.6, modified so West longitudes are negative and 0 is North
const gmst=greenwichMeanSiderealTime(jd_ut);
let localSiderealTime=(gmst+lon)%(2*Math.PI);

let H=(localSiderealTime - ra);
if(H<0){H+=2*Math.PI;}
if(H>Math.PI){H=H-2*Math.PI;}

let az = (Math.atan2(Math.sin(H), Math.cos(H)*Math.sin(lat) - Math.tan(dec)*Math.cos(lat)));
let a = (Math.asin(Math.sin(lat)*Math.sin(dec) + Math.cos(lat)*Math.cos(dec)*Math.cos(H)));
az-=Math.PI;

if(az<0){az+=2*Math.PI;}
return [az,a, localSiderealTime,H];
}

function greenwichMeanSiderealTime(jd){
//"Expressions for IAU 2000 precession quantities" N. Capitaine1,P.T.Wallace2, and J. Chapront
const t = ((jd - 2451545.0)) / 36525.0;

let gmst=this.earthRotationAngle(jd)+(0.014506 + 4612.156534*t + 1.3915817*t*t - 0.00000044 *t*t*t - 0.000029956*t*t*t*t - 0.0000000368*t*t*t*t*t)/60.0/60.0*Math.PI/180.0;  //eq 42
gmst%=2*Math.PI;
if(gmst<0) gmst+=2*Math.PI;

return gmst;
}

function earthRotationAngle(jd){
//IERS Technical Note No. 32

const t = jd- 2451545.0;
const f = jd%1.0;

let theta = 2*Math.PI * (f + 0.7790572732640 + 0.00273781191135448 * t); //eq 14
theta%=2*Math.PI;
if(theta<0)theta+=2*Math.PI;

return theta;

}


Example implementations are available: Sun Position RA/Dec to Alt/Az

• I just checked that declination formula against Horizons. It's quite accurate for the present era. Here's the Horizons data for this year Feb 14, 2022 at 18:19
• @greg-miller sorry I mentioned altitude instead of latitude and used along (wrong phrasing), my bad. Feb 15, 2022 at 5:55
• @PM2Ring an observer on the equator with a 0 degree decline of Sun will observe the 12 o'clock Sun at 0 degree N/S, also an observer at the tropic of cancer line will see the same or 23.45 degrees to south, the 12 o'clock Sun given on a specific day? Feb 15, 2022 at 9:02
• @Pavel The Sun has declination 0° on the equinoxes, for all observers. But its altitude angle (also known as the elevation) does depend on the observer's latitude. Please see en.wikipedia.org/wiki/Horizontal_coordinate_system Feb 15, 2022 at 12:08
• This answer shows how to calculate both the declination and altitude angle. It also computes the right ascension and azimuth angle. Feb 16, 2022 at 19:31