pidor/scripts/rscalc.awk

232 lines
5.5 KiB
Awk
Raw Normal View History

#!/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;
}