The following function sets the real sun position by setting time of day, day of year and current location (latitude) paramters:

Code:
#define DEGREE_IN_RADIANS 0.0174532925
#define EARTH_AXIAL_TILT 23.5			//earth axial tilt in degrees
#define EARTH_ORBITAL_PERIOD 365.25		//earth orbital period in days (length of a year)
#define EARTH_ROTATION_PERIOD 24.0		//earth rotation period in hours (length of a day)
#define LOCATION_TIME_SPRING_EQUINOX 81.0	//day number of the year when day and night cycle is exactly 12/12 hours
#define LOCATION_LATITUDE 48.0			//latitude of current location

void set_sun_position(var day_,var hour_,var minute_,var second_,var latitude_,var spring_equinox_)
{
	double B=(360.0/EARTH_ORBITAL_PERIOD) * (day_-spring_equinox_);

	//Equation of Time (EoT)
	double EoT=9.87*sin(2*B*DEGREE_IN_RADIANS) - 7.53*cos(B*DEGREE_IN_RADIANS) - 1.5*sin(B*DEGREE_IN_RADIANS);

	//Time Correction Factor (TC)
	double TC=EoT;

	//Local Time (LT)
	double LT=hour_ + (minute_/60.0) + (second_/3600.0);

	//Local Solar Time (LST)
	double LST=LT + TC/60.0;

	//Hour Angle (HRA)
	double HRA=15.0 * (LST - 12.0);

	//Number of days since the start of the year (d)
	double d=day_ + LT/24.0;

	//Declination (D)
	double D=EARTH_AXIAL_TILT * sin((360.0/EARTH_ORBITAL_PERIOD) * (d - spring_equinox_) * DEGREE_IN_RADIANS);

	//Sun Elevation (sun_angle.tilt)
	double sun_angle_tilt_sin=cos(HRA*DEGREE_IN_RADIANS)*cos(D*DEGREE_IN_RADIANS)*cos(LOCATION_LATITUDE*DEGREE_IN_RADIANS)+sin(D*DEGREE_IN_RADIANS)*sin(LOCATION_LATITUDE*DEGREE_IN_RADIANS);
	if(sun_angle_tilt_sin < -1){sun_angle_tilt_sin=-1;}
	if(sun_angle_tilt_sin > 1){sun_angle_tilt_sin=1;}		
	double sun_angle_tilt=asin(sun_angle_tilt_sin)/DEGREE_IN_RADIANS;

	//Sun Azimuth (sun_angle.pan)
	double sun_angle_pan_cos=(sin(D*DEGREE_IN_RADIANS)-sin(sun_angle_tilt*DEGREE_IN_RADIANS)*sin(LOCATION_LATITUDE*DEGREE_IN_RADIANS))/(cos(sun_angle_tilt*DEGREE_IN_RADIANS)*cos(LOCATION_LATITUDE*DEGREE_IN_RADIANS));
	if(sun_angle_pan_cos < -1){sun_angle_pan_cos=-1;}
	if(sun_angle_pan_cos > 1){sun_angle_pan_cos=1;}		
	double sun_angle_pan=acos(sun_angle_pan_cos)/DEGREE_IN_RADIANS;
	if(LST>12.0 && HRA>0.0){sun_angle_pan=360.0-sun_angle_pan;}
	sun_angle_pan=cycle(sun_angle_pan,0,360);
		
	sun_angle.pan=sun_angle_pan;
	sun_angle.tilt=sun_angle_tilt;
}



Here's an example how to call the function:
Code:
//set sun position to day 181 (june 30th) at time 14:30:00 (2:30 P.M.)
set_sun_position(181,14,30,0,LOCATION_LATITUDE,LOCATION_TIME_SPRING_EQUINOX);



The formulas of this function are taken from this websites:
http://pvcdrom.pveducation.org/SUNLIGHT/SUNPOS.HTM
http://pveducation.org/pvcdrom/properties-of-sunlight/sun-position-calculator
http://pveducation.org/pvcdrom/properties-of-sunlight/azimuth-angle
http://en.wikipedia.org/wiki/Solar_elevation_angle
http://en.wikipedia.org/wiki/Solar_azimuth_angle

Last edited by oliver2s; 12/24/13 21:14.