#define pi 3.141592653589793238462
struct time_struct
{
       int year;
       int month;
       int day;
       float hour;
       float min;
       float sec;
       float UTC;
};       

struct time_struct time;

struct position_struct
{
       float longitude;
       float latitude;
       float altitude;
};
struct position_struct position;

struct sun_position
{
       double zenith;
       double azimuth;
};
struct sun_position ENEA_OUTPUT;

double JD(struct time_struct time)
{
   double tg,day_fraction;
if (time.month<3)
{
   time.month = time.month + 12;
   time.year = time.year - 1;
}
day_fraction = (time.hour + time.min/60 + time.sec/3600)/24;

tg  = floor(365.25*(time.year - 2000)) + floor(30.6001*(time.month + 1))
    + time.day + time.UTC/24 + day_fraction - 1158.5;
   return tg;
}

struct sun_position ENEA(struct time_struct time, struct position_struct position)
{
//double output[2];
double delta_tt,theta_,fii,P,T,t,Sigma_,Ly,Lh,Lp,t2,Lm,JD2010;
double delta_psi, epsilon_, gamma_, alfa_, delta_, h, delta_alfa, alfa_t;
double delta_t, ht, cht, sht, e0, delta_e, z, A,L;   
JD2010 = JD(time) - 2556.5;

delta_tt = 66.18 - 0.0001033*cos(JD2010*0.0172) + 0.001654*sin(JD2010*0.0172);
theta_ = position.longitude*pi/180;
fii = position.latitude*pi/180;
P = 1;
T = 12;
t = JD(time) + delta_tt/86400;

// Calculation of the heliocentric longitude of the earth
Sigma_ = 1.72019e-2*t - 0.0563;


// Linear Term
Ly = 1.74094 + 1.720276868e-2*t + 3.34118e-2*sin(Sigma_) + 3.488e-4*sin(2*Sigma_);

// Moon Pertubation
Lm = 3.13e-5*sin(0.2127730*t - 0.585);
// Harmonic correction
Lh =  1.26e-5*sin(4.243e-3*t + 1.46);
Lh += 2.35e-5*sin(1.0727e-2*t + 0.72);
Lh += 2.76e-5*sin(1.5799e-2*t + 2.35);
Lh += 2.75e-5*sin(2.1551e-2*t - 1.98 );
Lh += 1.26e-5*sin(3.1490e-2*t - 0.80);
// Polynomial Correction

t2 = 0.001*t;
Lp = (((-2.3079e-7*t2 + 3.7976e-6)*t2 - 2.0458e-5)*t2 + 3.976e-5)*t2*t2;

L = Ly + Lm + Lh + Lp;

// Correction to geocentric longitude due to nutaion

delta_psi = 8.33e-5*sin(9.252e-4*t - 1.173);

// Earth axis inclination

epsilon_ = -6.21e-9*t + 0.409086 + 4.46e-5*sin(9.252e-4*t + 0.397);

// Geocentric global solar coordinates

// Geocentric solar longitude

gamma_ = L + pi + delta_psi - 9.932e-5;

// Geocentric right ascension

alfa_ = atan2(sin(gamma_)*cos(epsilon_),cos(gamma_));

// Geocentric Declination

delta_ = asin(sin(epsilon_)*sin(gamma_));

// Local hour angle of the sun

h = 6.30038809903*JD(time) + 4.8824623 + 0.9174*delta_psi +  theta_ - alfa_;

// Parallax correction to right ascension

delta_alfa = -4.26e-5*cos(fii)*cos(h);

// Topocentric sun coordinates

// Right Ascension

alfa_t = alfa_ + delta_alfa;

// Declination

delta_t = delta_ - 4.26e-5*(sin(fii) - delta_*cos(fii));

// Hour angle

ht = h - delta_alfa;
// ht = set_to_range(ht, 0, 2*pi);


// Approximate cosine and sine of ht
cht = cos(h) + delta_alfa*sin(h);
sht = sin(h) + delta_alfa*cos(h);

// Solar Elevation Angle without refraction correction

e0 = asin(sin(fii)*sin(delta_t) + cos(fii)*cos(delta_t)*cht);



// Atmospheric refraction correction to the solar elevation

delta_e = 0.084217*P/((273 + T)*tan(e0 + 0.0031376/(e0 + 0.089186)));

// Local Zenith
z = pi/2 - e0 - delta_e;
// Local Azimuth
A = atan2(sht,cht*sin(fii) - tan(delta_t)*cos(fii)) + pi;
ENEA_OUTPUT.zenith = 90 - z*180/pi;
ENEA_OUTPUT.azimuth = A*180/pi;
return ENEA_OUTPUT;
}