fn main
[
        if len(args)<5 then [
                write(h[2], nl+
                            "Usage: ajla sun.ajla longitude latitude timezone month day"+nl+
                            nl+
                            "Calculates sunset and sunrise (civil dawn and dusk) and solar noon in local time for given place, timezone and date."+nl+
                            nl+
                            "longitude in degrees.decimal_degree_fraction east  from Greenwich +/- east/west"  +nl+
                            "latitude  in degrees.decimal_degree_fraction north from equator   +/- north/south"+nl+
                            "timezone  in hours.decimal_hour_fraction     east  from Greenwich +/- east/west"+nl+
                            "month  1-12 integer"+nl+
                            "day    1-31 integer"+nl+
                            "Uses an imprecise criticised algorithm from Wikipedia https://en.wikipedia.org/wiki/Sunset_equation, see talk page."+nl+
                            "Wrong month or day number will throw an exception."+nl+
                            "In some places on some dates the sun doesn't set or rise."+nl+
                            nl);
		exit(1);
		return;
	]

	var M_PI :=instance_real_number_real64.pi;
	var lat_deg:real64:=instance_real_number_real64.from_bytes(args[0]);
	var lon_deg:real64:=instance_real_number_real64.from_bytes(args[1]);
	var tz_hours_eastward:=instance_real_number_real64.from_bytes(args[2]);
	var month:int:=ston(args[3]);
	var day:int:=ston(args[4]);

	{Incl. correction for aberration:
	16 min solar radius, 35 min horizon atmospheric refraction,
	6 deg or 360 from horizon to CENTER for civil twilight
	Sunset (top of Sun touches horizon):     -16-35
	Civil dusk (centre of Sun 6 deg below horizon): -35-360
	https://en.wikipedia.org/wiki/Sunset#/media/File:Twilight_subcategories.svg
	http://jgiesen.de/refract/img/refract.gif
	Navy algorithm (not used, used Wikipedia): https://aa.usno.navy.mil/faq/sun_approx
	}
	var solar_center_altitude_min:real64:= -35-360;

				{  31  28 31 30 31 30 31 31 30 31 30}
	var month_offsets:=list(int).[0, 1, -1, 0, 0, 1, 1, 2, 3, 3, 4, 4];
	var month_begins_days_since_jan1:int:=(30*(month-1))+month_offsets[month-1];
	var days_since_jan1:int:=month_begins_days_since_jan1+day-1;
	var lon_offset_h:=(0.+tz_hours_eastward)-(lon_deg/15.0);
	var Jstar2:real64:=0.0009+days_since_jan1;
	var M_deg:real64:=fmod(357.5291 + 0.98560028*Jstar2,360);
	var M_rad:real64:=M_deg*M_PI/180.0;
	var C:real64:=1.9148*sin(M_rad)+0.0200*sin(2.0*M_rad)+0.0003*sin(3.0*M_rad);
	var lambda_deg:real64:=fmod((M_deg+102.9372+C+180.0),360);
	var lambda_rad:real64:=lambda_deg*M_PI/180.0;
	var Jtransit_JD2:real64:= Jstar2 + (0.0053*sin(M_rad)) - (0.0069*sin(2.0*lambda_rad));
	var Jtransit_h:real64:=fmod(Jtransit_JD2,1)*24;
	var ax:=1;

	if Jtransit_h>=12.0 then Jtransit_h:=Jtransit_h-24.0;

	var decl_rad:real64:=asin(sin(lambda_rad)*sin(23.45*M_PI/180.0));
	var decl_deg:real64:=decl_rad*180.0/M_PI;
	var lat_rad:real64:=lat_deg/180.0*M_PI;
	var solar_center_altitude_deg:real64:=solar_center_altitude_min/60.0;
	var solar_center_altitude_rad:real64:=solar_center_altitude_deg*M_PI/180.0;
	var houra_rad:real64:=acos((sin(solar_center_altitude_rad)-sin(lat_rad)*sin(decl_rad))
					/(cos(lat_rad)*cos(decl_rad)));
	var houra_h:real64:=(houra_rad*12.0)/M_PI;
	var noon_h:real64:=12.0+Jtransit_h+lon_offset_h;
	var sunrise_h:real64:=noon_h-houra_h;
	var sunset_h:real64:=noon_h+houra_h;

	{ Add 1/2 min for rounding: }
	sunrise_h+=0.5/60.0;
	noon_h+=0.5/60.0;
	sunset_h+=0.5/60.0;
	var noon_m:=60.0*fract(noon_h);
	var sunrise_m:=60.0*fract(sunrise_h);
	var sunset_m:=60.0*fract(sunset_h);
	var sunrise_h_int:int:=floor(sunrise_h);
	var noon_h_int:int:=floor(noon_h);
	var sunset_h_int:int:=floor(sunset_h);
	var sunrise_m_int:int:=floor(sunrise_m);
	var noon_m_int:int:=floor(noon_m);
	var sunset_m_int:int:=floor(sunset_m);
	write(h[1], "Observer latitude is "+ ntos(lat_deg)+" deg from equator (+/- is north/south) and longitude "+ntos(lon_deg)
		+" deg east from Greenwich (+/- is east/west)."+nl);
	write(h[1], "Time zone is "+ ntos(tz_hours_eastward)+" hours from Greenwich (+/- is east/west)."+nl);
	write(h[1], "Month "+ ntos(month)+", day "+ntos(day) +nl);
	write(h[1], "Month "+ntos(month)+" begins "+ntos(month_begins_days_since_jan1)+" days since noon of Jan 1 (noon Jan is 0)."+nl);
	write(h[1], "Month "+ ntos(month)+", day "+ntos(day) +" is "+ntos(days_since_jan1)+" days since noon Jan 1 (noon Jan 1 is 0)."+nl);
	write(h[1], "Mean anomaly of the Sun is "+ ntos(M_deg)+" degrees since perihelion of Earth around January 3, counted increasing with time."+nl);
	write(h[1], "Equation of the center of the the Sun is "+ ntos(C)+" degrees."+nl);
	write(h[1], "Ecliptic longitude of the Sun is "+ ntos(lambda_deg)+" degrees, 0 occurs at spring equinox."+nl);
	write(h[1], "Solar noon (transit) is at "+ ntos(Jtransit_h)+" decimal hours relative to 12:00."+nl);
	write(h[1], "Declination of the Sun is "+ ntos(decl_deg)+" degrees north from equator."+nl);
	write(h[1], "Longitude offset in hours "+ntos(lon_offset_h) +nl);
	if not is_exception houra_h then [
		write(h[1], "Hour angle for civil dusk/dawn is +/- "+ ntos(houra_h)+" hours from the solar noon"+nl);
		write(h[1], "Civil dawn is at "+ ntos(sunrise_h)+" decimal hours ."+nl);
		write(h[1], "Solar noon is at "+ ntos(noon_h)+" decimal hours."+nl);
		write(h[1], "Civil dusk is at "+ ntos(sunset_h)+" decimal hours ."+nl);
		write(h[1], "Civil dawn is at "+ ntos(sunrise_h_int)+":"+list_left_pad(ntos(sunrise_m_int),2,'0')+" o'clock ."+nl);
		write(h[1], "Solar noon is at "+ ntos(noon_h_int)+":"+list_left_pad(ntos(noon_m_int),2,'0')+" o'clock ."+nl);
		write(h[1], "Civil dusk is at "+ ntos(sunset_h_int)+":"+list_left_pad(ntos(sunset_m_int),2,'0')+" o'clock ."+nl);
	] else [
		write(h[1], "Solar noon is at "+ ntos(noon_h)+" decimal hours."+nl);
		write(h[1], "Solar noon is at "+ ntos(noon_h_int)+":"+list_left_pad(ntos(noon_m_int),2,'0')+" o'clock ."+nl);
		write(h[1], "Civil dusk and dawn don't occur."+nl);
	]
]