Jump to content
Bohemia Interactive Forums
Sign in to follow this  
Followers 0
das attorney

Help needed: Calculating solar azimuth

By das attorney, in ARMA 2 & OA : MISSIONS - Editing & Scripting

Recommended Posts

das attorney    857

das attorney
Posted

I would like to calculate the solar azimuth angle for a mod I'm working on.

________________________________________________

________________________________________________

azelzen.gif

________________________________________________

________________________________________________

From this post by tpw, we can calculate the solar elevation angle (originally worked out by Carl Gustaffa):

_lat = -1 * getNumber(configFile >> "CfgWorlds" >> worldName >> "latitude"); //Arma latitude is negated for some odd reason.
_day = 360 * (dateToNumber date); //Convert current day to 360 for trigonometric calculations.
_hour = (daytime / 24) * 360; //Convert current hours to 360 for trigonometric calculations.
sunangle = ((12 * cos(_day) - 78) * cos(_lat) * cos(_hour)) - (24 * sin(_lat) * cos(_day));

I've tried following information on a few websites to calculate the solar azimuth angle, but I just can't get it to work. The following is code I've tried with the websites it's copied from with varying degrees of failure. Occasionally, I get the following value for the solar azimuth returned: -1#INT(?). If anyone can help in any way, I'd be very grateful.

NB: I put a lot of sleeps in to smooth out what's going on (while testing).

Method 1 (the script runs but returns wrong values including -1#INT whatever that is)

http://en.wikipedia.org/wiki/Solar_azimuth_angle

private ["_declin","_v1","_v2","_v3","_v4","_v5","_cos_azi","_azi","_lat","_day","_hour","_solar_elev_angle"];

// player switchCamera "GROUP";

while {true} do
{
_lat = -1 * getNumber(configFile >> "CfgWorlds" >> worldName >> "latitude");
_day = 360 * (dateToNumber date);
_day_2 = (360/365.25) * ( 365 * dateToNumber date);
_hour = (daytime / 24) * 360;
_solar_elev_angle = ((12 * cos(_day) - 78) * cos(_lat) * cos(_hour)) - (24 * sin(_lat) * cos(_day));

if (_solar_elev_angle > 0) then
{		
	_declin =  0.396372 - (22.91327 * cos(_day_2)) + (4.02543 * sin(_day_2)) - (0.387205 * cos(2 * _day_2)) + (0.051967 * sin(2 * _day_2)) - (0.154527 * cos(3 * _day_2)) + (0.084798 * sin(3 * _day_2));
	sleep 0.05;
	_v1 = sin _declin;
	sleep 0.05;
	_v2 = sin _solar_elev_angle;
	sleep 0.05;
	_v3 = sin _lat;
	sleep 0.05;
	_v4 = cos _solar_elev_angle;
	sleep 0.05;
	_v5 = cos _lat;
	sleep 0.05;
	_v6 = getDir player;
	sleep 0.05;

	_f1 = _v1 - _v2;
	sleep 0.05;
	_f2 = _f1 * _v3;
	sleep 0.05;
	_f3 = _v4 * _v5;
	sleep 0.05;
	_cos_azi = _f2 / _f3;
	sleep 0.05;

	_azi = acos _cos_azi;

	hint format ["azi: %1\nsin_Dec: %2\nsin_SEA: %3\nsin_lat: %4\ncos_SEA: %5\ncos_lat: %6\nPlayer dir: %7", _azi, _v1, _v2, _v3, _v4, _v5, _v6];
	hint format ["azi: %1\nDec: %2\nSEA: %3\nlat: %4", _azi, _declin, _solar_elev_angle, _lat];

	sleep 0.5;

};
skipTime 0.2;
sleep 0.1;
};

Method 2 (the script runs but returns wrong values)

http://answers.google.com/answers/threadview/id/782886.html

private ["_sha","_cos_sza","_dec","_tc","_sza","_cos_azi","_azi","_lat","_long","_day","_dayz","_hour","_solar_elev_angle"];
while {true} do
{
_lat = -1 * getNumber(configFile >> "CfgWorlds" >> worldName >> "latitude");
_long = getNumber(configFile >> "CfgWorlds" >> worldName >> "longitude");
_dayz = 360 * (dateToNumber date);
_day = (360/365.25) * ( 365 * dateToNumber date);
_hourz = (daytime / 24) * 360;
_hour = daytime;
_solar_elev_angle = ((12 * cos(_dayz) - 78) * cos(_lat) * cos(_hourz)) - (24 * sin(_lat) * cos(_dayz));

if (_solar_elev_angle > 0) then
{
	_dec =  0.396372 - (22.91327 * cos(_day)) + (4.02543 * sin(_day)) - (0.387205 * cos(2 * _day)) + (0.051967 * sin(2 * _day)) - (0.154527 * cos(3 * _day)) + (0.084798 * sin(3 * _day));

	sleep 0.05;

	_tc = 0.004297 + (0.107029 * cos(_day)) - (1.837877 * sin(_day)) - (0.837378 * cos(2*_day)) - (2.340475 * sin(2*_day));

	sleep 0.05;

	_sha = ((_hour - 12) * 15) + _long + _tc;

	if (_sha > 180) then {_sha = _sha -360};
	if (_sha < -180) then {_sha = _sha + 360};

       _cos_sza = (sin(_lat)) * (sin(_dec)) + (cos(_lat)) * (cos(_dec)) * (cos(_sha));

	sleep 0.05;

	if (_cos_sza > 1) then {_cos_sza = 1};
	if (_cos_sza < -1) then {_cos_sza = -1};		

	_sza = acos _cos_sza;

	sleep 0.05;

       _cos_azi = ((sin(_dec)) - (sin(_lat)) * (cos(_sza))) / ((cos(_lat)) * (sin(_sza)));

	sleep 0.05;

	_azi = acos _cos_azi;

	hint format ["azi: %1\nDec: %2\nTC: %3\nSHA: %4\nSZA: %5", _azi, _dec, _tc, _sha, _sza];

	sleep 0.5;
};

skipTime 0.2;
sleep 0.1;
};

Method 3 (the script runs but returns wrong values including -1#INT)

http://www.jgiesen.de/astro/suncalc/calculations.htm

private ["_sha","_cos_phi","_phi","_azi","_lat","_dayz","_solar_elev_angle","_year","_eqtime","_declin","_offset","_time_offset","_tst","_cos180_azi","_180_azi","_lon","_hourz"];
while {true} do
{
_lat = getNumber(configFile >> "CfgWorlds" >> worldName >> "latitude");
_lon = getNumber(configFile >> "CfgWorlds" >> worldName >> "longitude");
_dayz = 360 * (dateToNumber date);
_hourz = (daytime / 24) * 360;
_solar_elev_angle = ((12 * cos(_dayz) - 78) * cos(_lat) * cos(_hourz)) - (24 * sin(_lat) * cos(_dayz));

if (_solar_elev_angle > 0) then
{		
	// fractional year in rads
       _year = (2 * pi / 365) * (((floor _dayz) - 1) + ((floor daytime) - 12) / 24);

	// equation of time
       _eqtime = 229.18 * (0.000075 + (0.001868 * (cos _year)) - (0.032077 * (sin _year)) - (0.014615 * (cos(2 *_year))) - (0.040849 * (sin( 2 * _year))));

	_declin = 0.006918 - (0.399912 * (cos _year)) + (0.070257 * (sin _year)) - (0.006758 * (cos(2 * _year))) + (0.000907 * (sin(2 * _year))) - (0.002697 * (cos(3 * _year))) + (0.00148 * (sin(3 *_year)));

	_offset = _lon / 15;

	_time_offset = (_eqtime - 4) * _lon + 60 * _offset;

	// true solar time in mins
       _tst = (round(daytime * 60)) + _time_offset;

	// solar hour angle in deg
	_sha = (_tst / 4) - 180;

	// solar zenith angle (phi)

       _cos_phi = ((sin _lat) * (sin _declin) + (cos _lat) * (cos _declin) * (cos _sha));

	if (_cos_phi > 1) then {_cos_phi = 1};
	if (_cos_phi < -1) then {_cos_phi = -1};

	_phi = acos _cos_phi;

       _cos180_azi = -1 * ((((sin _lat) * (cos _phi)) - (sin _declin)) / (cos _lat) * (sin _phi));

	_180_azi = acos _cos180_azi;

	_azi = _180_azi + 180;

	hint format ["azi: %1", _azi];

	sleep 0.5;

};

skipTime 0.2;
sleep 0.1;
};

Share this post

Link to post
Share on other sites

das attorney    857

das attorney
Posted (edited)

Been working on this a bit and found a script that 'sort of' works.

The suns position is mainly correct but can be off by 30 deg at sunset and sunrise. Sometimes I'm still getting an invalid number returned.

I didn't really expect any replies but still glad for help if available.

Source:

http://pveducation.org/pvcdrom/properties-of-sunlight/suns-position

Script:

private ["_cos_azi","_azi","_azi1","_lat","_day","_sea_day","_sea_hour","_sea","_lstm","_b","_eot","_tc","_lst","_hra","_dec","_max_elev_at_solar_noon","_lon","_time_zone","_pos","_x_pos","_diff","_fraction","_map_lat","_center_pos","_half_width","_width"];

_lat = -1 * getNumber(configFile >> "CfgWorlds" >> worldName >> "latitude");
_lon = getNumber(configFile >> "CfgWorlds" >> worldName >> "longitude");
_time_zone = 2; // Chernarus

while {true} do
{
_day = round floor( 365 * dateToNumber date);
_sea_day = 360 * (dateToNumber date);
_sea_hour = (daytime / 24) * 360;
_sea = ((12 * cos(_sea_day) - 78) * cos(_lat) * cos(_sea_hour)) - (24 * sin(_lat) * cos(_sea_day));

// _sea = ASin [ Sin _dec * Sin _lat + Cos _dec * Cos _hra * Cos _lat];
if (_sea > 0) then
{		
	_lstm = 15 * _time_zone;
	// _lstm = _lon; // works as well..

	_b = (360 / 365) * (_day - 81);

	_eot = 9.87 * sin(2 * _b) - 7.53 * cos(_b) - 1.5 * sin(_b);

	_tc = 4 * (_lon - _lstm) + _eot;

	_lst = daytime + (_tc / 60);

	_hra = 15 * (_lst - 12);

	_dec = 23.45 * sin(_b);

	// northern hemishere !!  // _max_elev_at_solar_noon = 90 + _lat - _dec; for southern hemi

	_max_elev_at_solar_noon = 90 - _lat + _dec;

       _cos_azi = (((sin(_dec)) * cos(_lat)) - ((cos(_dec)) * sin(_lat) * cos(_hra))) / cos _max_elev_at_solar_noon;

	if (_lst < 12 or _hra < 0) then // morning
	{
		_azi1 = acos _cos_azi;
	} else // afternoon
	{
		_azi1 = 360 - (acos _cos_azi);
	};

	diag_log format ["time: %1 azi: %2", daytime, _azi1];

	hint format ["azi: %1\nDec: %2\nHRA: %3\nLST: %4\nLSTM: %5", _azi1, _dec, _hra, _lst, _lstm];

	// player setDir _azi1;

	sleep 0.5;

};

skipTime 0.2;
sleep 0.1;
};

Please see following for rpt readouts of the azimuth vs time:

// Chernarus May 6, 2013

// "time: 5.53457 azi: 56.1218"

// "time: 5.73477 azi: 61.0567"

// "time: 5.93496 azi: 65.7796"

// "time: 6.13516 azi: 70.3337"

// "time: 6.33535 azi: 74.7502"

// "time: 6.53555 azi: 79.0524"

// "time: 6.73574 azi: 83.2587"

// "time: 6.93594 azi: 87.3829"

// "time: 7.13614 azi: 91.4369"

// "time: 7.33633 azi: 95.4297"

// "time: 7.53652 azi: 99.3692"

// "time: 7.73672 azi: 103.262"

// "time: 7.93692 azi: 107.113"

// "time: 8.13711 azi: 110.928"

// "time: 8.33731 azi: 114.711"

// "time: 8.53751 azi: 118.464"

// "time: 8.7377 azi: 122.191"

// "time: 8.93789 azi: 125.896"

// "time: 9.13809 azi: 129.579"

// "time: 9.33828 azi: 133.245"

// "time: 9.53848 azi: 136.894"

// "time: 9.73867 azi: 140.528"

// "time: 9.93886 azi: 144.149"

// "time: 10.1391 azi: 147.76"

// "time: 10.3393 azi: 151.36"

// "time: 10.5395 azi: 154.951"

// "time: 10.7397 azi: 158.535"

// "time: 10.9399 azi: 162.113"

// "time: 11.1401 azi: 165.685"

// "time: 11.3403 azi: 169.253"

// "time: 11.5405 azi: 172.818"

// "time: 11.7407 azi: 176.381"

// "time: 11.9409 azi: 179.937"

// "time: 12.1411 azi: 183.505"

// "time: 12.3413 azi: 187.068"

// "time: 12.5415 azi: 190.633"

// "time: 12.7417 azi: 194.201"

// "time: 12.9419 azi: 197.774"

// "time: 13.1421 azi: 201.351"

// "time: 13.3424 azi: 204.935"

// "time: 13.5426 azi: 208.526"

// "time: 13.7428 azi: 212.125"

// "time: 13.943 azi: 215.735"

// "time: 14.1432 azi: 219.356"

// "time: 14.3434 azi: 222.99"

// "time: 14.5435 azi: 226.638"

// "time: 14.7437 azi: 230.303"

// "time: 14.9439 azi: 233.986"

// "time: 15.1441 azi: 237.69"

// "time: 15.3443 azi: 241.416"

// "time: 15.5445 azi: 245.169"

// "time: 15.7447 azi: 248.95"

// "time: 15.9449 azi: 252.764"

// "time: 16.1451 azi: 256.614"

// "time: 16.3453 azi: 260.505"

// "time: 16.5455 azi: 264.443"

// "time: 16.7457 azi: 268.434"

// "time: 16.9459 azi: 272.486"

// "time: 17.1461 azi: 276.608"

// "time: 17.3463 azi: 280.811"

// "time: 17.5465 azi: 285.11"

// "time: 17.7467 azi: 289.523"

// "time: 17.9469 azi: 294.072"

// "time: 18.1471 azi: 298.789"

// "time: 18.3473 azi: 303.716"

// "time: 18.5474 azi: 308.915"

// Chernarus June 6, 2013

// "time: 5.1362 azi: -1.#IND"

// "time: 5.33639 azi: 6.28051"

// "time: 5.53658 azi: 25.1475"

// "time: 5.73678 azi: 35.3234"

// "time: 5.93698 azi: 43.4232"

// "time: 6.13717 azi: 50.4574"

// "time: 6.33737 azi: 56.8243"

// "time: 6.53757 azi: 62.728"

// "time: 6.73776 azi: 68.2895"

// "time: 6.93796 azi: 73.5873"

// "time: 7.13815 azi: 78.6754"

// "time: 7.33835 azi: 83.5928"

// "time: 7.53854 azi: 88.3686"

// "time: 7.73874 azi: 93.0257"

// "time: 7.93893 azi: 97.5817"

// "time: 8.13913 azi: 102.051"

// "time: 8.33932 azi: 106.445"

// "time: 8.53952 azi: 110.774"

// "time: 8.73971 azi: 115.046"

// "time: 8.93991 azi: 119.267"

// "time: 9.14013 azi: 123.445"

// "time: 9.34035 azi: 127.585"

// "time: 9.54056 azi: 131.689"

// "time: 9.74077 azi: 135.764"

// "time: 9.94098 azi: 139.812"

// "time: 10.1412 azi: 143.837"

// "time: 10.3414 azi: 147.841"

// "time: 10.5416 azi: 151.828"

// "time: 10.7418 azi: 155.8"

// "time: 10.942 azi: 159.759"

// "time: 11.1422 azi: 163.708"

// "time: 11.3424 azi: 167.649"

// "time: 11.5425 azi: 171.583"

// "time: 11.7427 azi: 175.514"

// "time: 11.9429 azi: 179.442"

// "time: 12.1431 azi: 183.369"

// "time: 12.3433 azi: 187.299"

// "time: 12.5435 azi: 191.232"

// "time: 12.7436 azi: 195.171"

// "time: 12.9438 azi: 199.117"

// "time: 13.144 azi: 203.073"

// "time: 13.3442 azi: 207.041"

// "time: 13.5444 azi: 211.024"

// "time: 13.7446 azi: 215.023"

// "time: 13.9449 azi: 219.042"

// "time: 14.1451 azi: 223.083"

// "time: 14.3453 azi: 227.15"

// "time: 14.5455 azi: 231.246"

// "time: 14.7457 azi: 235.374"

// "time: 14.9459 azi: 239.54"

// "time: 15.1461 azi: 243.749"

// "time: 15.3463 azi: 248.005"

// "time: 15.5465 azi: 252.317"

// "time: 15.7467 azi: 256.692"

// "time: 15.9469 azi: 261.139"

// "time: 16.1471 azi: 265.669"

// "time: 16.3473 azi: 270.295"

// "time: 16.5475 azi: 275.036"

// "time: 16.7477 azi: 279.91"

// "time: 16.9478 azi: 284.946"

// "time: 17.148 azi: 290.18"

// "time: 17.3482 azi: 295.659"

// "time: 17.5484 azi: 301.456"

// "time: 17.7486 azi: 307.675"

// "time: 17.9488 azi: 314.491"

// "time: 18.149 azi: 322.226"

// "time: 18.3492 azi: 331.625"

// "time: 18.5494 azi: 345.653"

// "time: 18.7496 azi: -1.#IND"

// Chernarus July 6, 2013

// "time: 5.13784 azi: -1.#IND"

// "time: 5.33804 azi: -1.#IND"

// "time: 5.53823 azi: 14.9756"

// "time: 5.73843 azi: 28.8203"

// "time: 5.93862 azi: 38.203"

// "time: 6.13882 azi: 45.9448"

// "time: 6.33901 azi: 52.774"

// "time: 6.53921 azi: 59.0089"

// "time: 6.7394 azi: 64.8221"

// "time: 6.9396 azi: 70.319"

// "time: 7.1398 azi: 75.5695"

// "time: 7.33999 azi: 80.6224"

// "time: 7.54018 azi: 85.5137"

// "time: 7.74038 azi: 90.2702"

// "time: 7.94058 azi: 94.9133"

// "time: 8.14077 azi: 99.4593"

// "time: 8.341 azi: 103.923"

// "time: 8.54122 azi: 108.314"

// "time: 8.74143 azi: 112.641"

// "time: 8.94165 azi: 116.914"

// "time: 9.14186 azi: 121.138"

// "time: 9.34207 azi: 125.32"

// "time: 9.54228 azi: 129.464"

// "time: 9.74248 azi: 133.575"

// "time: 9.94268 azi: 137.657"

// "time: 10.1429 azi: 141.713"

// "time: 10.3431 azi: 145.747"

// "time: 10.5433 azi: 149.762"

// "time: 10.7435 azi: 153.76"

// "time: 10.9437 azi: 157.743"

// "time: 11.1438 azi: 161.715"

// "time: 11.344 azi: 165.678"

// "time: 11.5442 azi: 169.633"

// "time: 11.7444 azi: 173.584"

// "time: 11.9446 azi: 177.53"

// "time: 12.1448 azi: 181.476"

// "time: 12.3449 azi: 185.422"

// "time: 12.5451 azi: 189.371"

// "time: 12.7453 azi: 193.325"

// "time: 12.9455 azi: 197.286"

// "time: 13.1457 azi: 201.255"

// "time: 13.3459 azi: 205.236"

// "time: 13.5461 azi: 209.23"

// "time: 13.7463 azi: 213.24"

// "time: 13.9465 azi: 217.269"

// "time: 14.1467 azi: 221.32"

// "time: 14.3469 azi: 225.395"

// "time: 14.5471 azi: 229.498"

// "time: 14.7473 azi: 233.633"

// "time: 14.9475 azi: 237.805"

// "time: 15.1477 azi: 242.018"

// "time: 15.3479 azi: 246.278"

// "time: 15.5481 azi: 250.592"

// "time: 15.7483 azi: 254.965"

// "time: 15.9485 azi: 259.409"

// "time: 16.1487 azi: 263.932"

// "time: 16.3489 azi: 268.55"

// "time: 16.5491 azi: 273.276"

// "time: 16.7493 azi: 278.131"

// "time: 16.9495 azi: 283.14"

// "time: 17.1497 azi: 288.337"

// "time: 17.3499 azi: 293.766"

// "time: 17.5501 azi: 299.491"

// "time: 17.7503 azi: 305.607"

// "time: 17.9505 azi: 312.264"

// "time: 18.1507 azi: 319.732"

// "time: 18.3509 azi: 328.584"

// "time: 18.5511 azi: 340.636"

// "time: 18.7513 azi: -1.#IND"

// "time: 18.9514 azi: -1.#IND"

// Chernarus August 6, 2013

// "time: 5.13498 azi: 37.7085"

// "time: 5.33518 azi: 44.122"

// "time: 5.53538 azi: 49.9238"

// "time: 5.73557 azi: 55.307"

// "time: 5.93576 azi: 60.3841"

// "time: 6.13596 azi: 65.2273"

// "time: 6.33615 azi: 69.8852"

// "time: 6.53635 azi: 74.3934"

// "time: 6.73655 azi: 78.7778"

// "time: 6.93675 azi: 83.0585"

// "time: 7.13694 azi: 87.2512"

// "time: 7.33714 azi: 91.3682"

// "time: 7.53733 azi: 95.4199"

// "time: 7.73753 azi: 99.4146"

// "time: 7.93773 azi: 103.36"

// "time: 8.13792 azi: 107.26"

// "time: 8.33812 azi: 111.122"

// "time: 8.53831 azi: 114.95"

// "time: 8.73851 azi: 118.747"

// "time: 8.93873 azi: 122.517"

// "time: 9.13895 azi: 126.262"

// "time: 9.33917 azi: 129.986"

// "time: 9.53938 azi: 133.69"

// "time: 9.73959 azi: 137.377"

// "time: 9.93979 azi: 141.048"

// "time: 10.14 azi: 144.705"

// "time: 10.3402 azi: 148.351"

// "time: 10.5404 azi: 151.986"

// "time: 10.7406 azi: 155.612"

// "time: 10.9408 azi: 159.23"

// "time: 11.141 azi: 162.842"

// "time: 11.3412 azi: 166.448"

// "time: 11.5413 azi: 170.05"

// "time: 11.7415 azi: 173.65"

// "time: 11.9417 azi: 177.247"

// "time: 12.1419 azi: 180.843"

// "time: 12.3421 azi: 184.44"

// "time: 12.5423 azi: 188.038"

// "time: 12.7424 azi: 191.638"

// "time: 12.9426 azi: 195.242"

// "time: 13.1428 azi: 198.85"

// "time: 13.343 azi: 202.465"

// "time: 13.5432 azi: 206.087"

// "time: 13.7434 azi: 209.717"

// "time: 13.9436 azi: 213.357"

// "time: 14.1438 azi: 217.008"

// "time: 14.344 azi: 220.672"

// "time: 14.5442 azi: 224.35"

// "time: 14.7444 azi: 228.044"

// "time: 14.9446 azi: 231.757"

// "time: 15.1449 azi: 235.491"

// "time: 15.3451 azi: 239.247"

// "time: 15.5453 azi: 243.029"

// "time: 15.7455 azi: 246.84"

// "time: 15.9457 azi: 250.683"

// "time: 16.1459 azi: 254.563"

// "time: 16.3461 azi: 258.484"

// "time: 16.5462 azi: 262.451"

// "time: 16.7464 azi: 266.472"

// "time: 16.9466 azi: 270.553"

// "time: 17.1468 azi: 274.704"

// "time: 17.347 azi: 278.936"

// "time: 17.5472 azi: 283.263"

// "time: 17.7474 azi: 287.703"

// "time: 17.9476 azi: 292.278"

// "time: 18.1478 azi: 297.017"

// "time: 18.348 azi: 301.963"

// "time: 18.5482 azi: 307.173"

// "time: 18.7484 azi: 312.734"

Edited by Das Attorney

Share this post

Link to post
Share on other sites

panther42    40

panther42
Posted

Das Attorney, this stuff is way above me, but ruebe has a function in his library called fn_sun.sqf

Looks alot like some of the calculations you're doing here.

   Author:
   rübe

  Description:
   calculates sunrise/sunset on a given date, 
   based on the algorithm from:

    Almanac for Computers, 1990
    published by Nautical Almanac Office
    United States Naval Observatory
    Washington, DC 20392

    (see williams.best.vwh.net/sunrise_sunset_algorithm.htm)

Might be of help...

Share this post

Link to post
Share on other sites

Please sign in to comment

You will be able to leave a comment after signing in



Sign In Now
Sign in to follow this  
Followers 0
Go To Topic Listing ARMA 2 & OA : MISSIONS - Editing & Scripting
×