232 lines
5.5 KiB
Awk
232 lines
5.5 KiB
Awk
|
#!/usr/bin/awk -f
|
||
|
# converted from C to awk by gunstick@syn2cat.lu 2014/10/02
|
||
|
# based on http://www.sci.fi/~benefon/rscalc_cpp.html
|
||
|
# ugly C code gives ugly awk code, don't blame me
|
||
|
# input is like this: 49.5916 6.1407 2
|
||
|
|
||
|
# gawk program calculating the sunrise and sunset for
|
||
|
# the current date and a fixed location(latitude,longitude)
|
||
|
# Note, twilight calculation gives insufficient accuracy of results
|
||
|
# Jarmo Lammi 1999 - 2001
|
||
|
# Last update July 21st, 2001
|
||
|
|
||
|
BEGIN {
|
||
|
|
||
|
pi = 3.14159;
|
||
|
degs = 0.0;
|
||
|
rads = 0.0;
|
||
|
|
||
|
L=g=daylen=0.0
|
||
|
SunDia = 0.53; # Sunradius degrees
|
||
|
|
||
|
AirRefr = 34.0 / 60.0; # athmospheric refraction degrees
|
||
|
}
|
||
|
# Get the days to J2000
|
||
|
# h is UT in decimal hours
|
||
|
# FNday only works between 1901 to 2099 - see Meeus chapter 7
|
||
|
|
||
|
function FNday (y, m, d, h) {
|
||
|
luku = int(-7 * (y + (m + 9)/12)/4 + 275*m/9 + d);
|
||
|
# type casting necessary on PC DOS and TClite to avoid overflow
|
||
|
luku+= y*367;
|
||
|
return luku - 730531.5 + h/24.0;
|
||
|
};
|
||
|
|
||
|
# the function below returns an angle in the range
|
||
|
# 0 to 2*pi
|
||
|
|
||
|
function FNrange ( x) {
|
||
|
b = 0.5*x / pi;
|
||
|
a = 2.0*pi * (b - int(b));
|
||
|
if (a < 0) a = 2.0*pi + a;
|
||
|
return a;
|
||
|
};
|
||
|
|
||
|
# Calculating the hourangle
|
||
|
#
|
||
|
function f0( lat, declin) {
|
||
|
fo=dfo=0.0;
|
||
|
# Correction: different sign at S HS
|
||
|
dfo = rads*(0.5*SunDia + AirRefr)
|
||
|
if (lat < 0.0) { dfo = -dfo };
|
||
|
fo = tan(declin + dfo) * tan(lat*rads);
|
||
|
if (fo>0.99999) { fo=1.0 }; # to avoid overflow //
|
||
|
fo = asin(fo) + pi/2.0;
|
||
|
return fo;
|
||
|
};
|
||
|
|
||
|
# Calculating the hourangle for twilight times
|
||
|
#
|
||
|
function f1(lat, declin) {
|
||
|
fi=df1=0.0;
|
||
|
# Correction: different sign at S HS
|
||
|
df1 = rads * 6.0; if (lat < 0.0) { df1 = -df1 };
|
||
|
fi = tan(declin + df1) * tan(lat*rads);
|
||
|
if (fi>0.99999) { fi=1.0 } ; # to avoid overflow //
|
||
|
fi = asin(fi) + pi/2.0;
|
||
|
return fi;
|
||
|
};
|
||
|
|
||
|
|
||
|
# Find the ecliptic longitude of the Sun
|
||
|
|
||
|
function FNsun (d) {
|
||
|
|
||
|
# mean longitude of the Sun
|
||
|
|
||
|
L = FNrange(280.461 * rads + .9856474 * rads * d);
|
||
|
|
||
|
# mean anomaly of the Sun
|
||
|
|
||
|
g = FNrange(357.528 * rads + .9856003 * rads * d);
|
||
|
|
||
|
# Ecliptic longitude of the Sun
|
||
|
|
||
|
return FNrange(L + 1.915 * rads * sin(g) + .02 * rads * sin(2 * g));
|
||
|
};
|
||
|
|
||
|
|
||
|
# Display decimal hours in hours and minutes
|
||
|
function showhrmn( dhr) {
|
||
|
hr=mn=0;
|
||
|
hr= int(dhr);
|
||
|
mn = int((dhr - hr)*60);
|
||
|
|
||
|
printf("%0d:%0d",hr,mn);
|
||
|
};
|
||
|
|
||
|
# awk misses some trigo functions
|
||
|
function asin(x) { return atan2(x, sqrt(1-x*x)) }
|
||
|
function tan(x) { return sin(x)/cos(x) }
|
||
|
{
|
||
|
|
||
|
y=m=day=h=latit=longit=0.0
|
||
|
inlat=inlon=intz=0
|
||
|
tzone=d=lambda=0
|
||
|
obliq=alpha=delta=LL=equation=ha=hb=twx=0
|
||
|
twam=altmax=noont=settm=riset=twpm=0
|
||
|
sekunnit=0;
|
||
|
tm =0;
|
||
|
|
||
|
degs = 180.0/pi;
|
||
|
rads = pi/180.0;
|
||
|
# get the date and time from the user
|
||
|
# read system date and extract the year
|
||
|
|
||
|
# First get time **/
|
||
|
sekunnit=systime();
|
||
|
|
||
|
# Next get localtime **/
|
||
|
|
||
|
#p=localtime(&sekunnit);
|
||
|
|
||
|
y = strftime("%Y",sekunnit) # p->tm_year;
|
||
|
# this is Y2K compliant method
|
||
|
#y+= 1900;
|
||
|
m= strftime("%m",sekunnit) #m = p->tm_mon + 1;
|
||
|
|
||
|
|
||
|
day=strftime("%d",sekunnit) # day = p->tm_mday;
|
||
|
|
||
|
h = 12;
|
||
|
|
||
|
printf("year %4d month %2d\n",y,m);
|
||
|
printf("Input latitude, longitude [and timezone]\n");
|
||
|
#scanf("%f", &inlat); scanf("%f", &inlon);
|
||
|
#scanf("%f", &intz);
|
||
|
inlat=$1 ; inlon=$2
|
||
|
if($3 != "") {
|
||
|
intz=$3
|
||
|
} else { # guess the timezone
|
||
|
intz=(systime()-mktime(strftime("%Y %m %d %H %M %S",systime(),1)))/3600
|
||
|
}
|
||
|
latit = inlat; longit = inlon;
|
||
|
tzone = intz;
|
||
|
|
||
|
# testing
|
||
|
# m=6; day=10;
|
||
|
|
||
|
|
||
|
d = FNday(y, m, day, h);
|
||
|
|
||
|
# Use FNsun to find the ecliptic longitude of the
|
||
|
# Sun
|
||
|
|
||
|
lambda = FNsun(d);
|
||
|
|
||
|
# Obliquity of the ecliptic
|
||
|
|
||
|
obliq = 23.439 * rads - .0000004 * rads * d;
|
||
|
|
||
|
# Find the RA and DEC of the Sun
|
||
|
|
||
|
alpha = atan2(cos(obliq) * sin(lambda), cos(lambda));
|
||
|
delta = asin(sin(obliq) * sin(lambda));
|
||
|
|
||
|
# Find the Equation of Time
|
||
|
# in minutes
|
||
|
# Correction suggested by David Smith
|
||
|
LL = L - alpha;
|
||
|
if (L < pi) { LL += 2.0*pi } ;
|
||
|
equation = 1440.0 * (1.0 - LL / pi/2.0);
|
||
|
ha = f0(latit,delta);
|
||
|
hb = f1(latit,delta);
|
||
|
twx = hb - ha; # length of twilight in radians
|
||
|
twx = 12.0*twx/pi; # length of twilight in hours
|
||
|
printf "ha= %.2f hb= %.2f \n",ha,hb ;
|
||
|
# Conversion of angle to hours and minutes //
|
||
|
daylen = degs*ha/7.5;
|
||
|
if (daylen<0.0001) {daylen = 0.0;}
|
||
|
# arctic winter //
|
||
|
|
||
|
riset = 12.0 - 12.0 * ha/pi + tzone - longit/15.0 + equation/60.0;
|
||
|
settm = 12.0 + 12.0 * ha/pi + tzone - longit/15.0 + equation/60.0;
|
||
|
noont = riset + 12.0 * ha/pi;
|
||
|
altmax = 90.0 + delta * degs - latit;
|
||
|
# Correction for S HS suggested by David Smith
|
||
|
# to express altitude as degrees from the N horizon
|
||
|
if (latit < delta * degs) altmax = 180.0 - altmax;
|
||
|
|
||
|
twam = riset - twx; # morning twilight begin
|
||
|
twpm = settm + twx; # evening twilight end
|
||
|
|
||
|
if (riset > 24.0) riset-= 24.0;
|
||
|
if (settm > 24.0) settm-= 24.0;
|
||
|
|
||
|
print "\n sunrise and set";
|
||
|
print "===============";
|
||
|
|
||
|
printf " year : %d \n",y ;
|
||
|
printf " month : %d \n",m ;
|
||
|
printf " day : %d \n\n",day ;
|
||
|
printf "Days since Y2K : %d \n",d;
|
||
|
|
||
|
printf "Latitude : %3.1f, longitude: %3.1f, timezone: %3.1f \n",latit,longit,tzone ;
|
||
|
printf "Declination : %.2f \n",delta * degs ;
|
||
|
printf "Daylength : " ; showhrmn(daylen); printf " hours \n" ;
|
||
|
print ""
|
||
|
printf "Civil twilight: " ;
|
||
|
showhrmn(twam);
|
||
|
print ""
|
||
|
printf "Sunrise : " ;
|
||
|
showhrmn(riset);
|
||
|
printf " %d",mktime(strftime("%Y %m %d 00 00 00"))+int(riset*3600)
|
||
|
print ""
|
||
|
|
||
|
printf "Sun altitude " ;
|
||
|
# Amendment by D. Smith
|
||
|
printf " %.2f degr",altmax ;
|
||
|
printf latit>=0.0 ? " South" : " North" ;
|
||
|
printf " at noontime " ; showhrmn(noont); ;
|
||
|
print ""
|
||
|
printf "Sunset : " ;
|
||
|
showhrmn(settm);
|
||
|
printf " %d",mktime(strftime("%Y %m %d 00 00 00"))+int(settm*3600)
|
||
|
print ""
|
||
|
printf "Civil twilight: " ;
|
||
|
showhrmn(twpm); print "" ;
|
||
|
print ""
|
||
|
|
||
|
#return 0;
|
||
|
}
|