#!/usr/bin/ksh

# ksh or binfmtc ? php not really C and uses float rather than double
# abs acos acosh asin asinh atan atan2 atanh cbrt copysign cos cosh erf
# erfc exp exp2 expm1 fabs fdim finite floor fma fmax fmod hypot ilogb
# int isinf isnan lgamma log log2 logb nearbyint nextafter nexttoward pow
# remainder rint round sin sinh sqrt tan tanh tgamma trunc

# http://www.bom.gov.au/climate/environ/housedesign/solar_access.shtml

# Perth Astronomical Observatory, Bickley, on (09) 9293 8109 or email
# http://www.perthobservatory.wa.gov.au/ mailto:perthobs@iinet.au
# mailto:Perth.Observatory@dec.wa.gov.au
# Ask nicely they may check/comment on calcs

#http://www.usc.edu/dept/architecture/mbs/tools/vrsolar/Help/solar_concepts.html
# Archit Jain VRSolar mailto:architj@usc.edu
# Archit Jain(1999) Interactive Teaching Tool on the Web: Site Solar Access 
# A. Jain (1997) VRSolar:
#  An Exploration in Web Based Interactive Architectural Teaching
# Only some Java like frags available

# Wikipedia: --
# http://en.wikipedia.org/wiki/Declination#Sun
# The formula neglects the effect of atmospheric refraction.
# The formula is based on geocentric latitude, which differs from
# geographic latitude by less than 12 minutes of arc.
# http://en.wikipedia.org/wiki/Solar_elevation_angle
# http://en.wikipedia.org/wiki/Solar_azimuth_angle
# Tidy up and structure with
# http://en.wikipedia.org/wiki/Trigonometric_identities

# http://www.moshier.net/ radiation calcs
# Redo AS 2627.1 with correct rather than close solution
# http://van-der-waals.pc.uni-koeln.de/quartic/quartic.html

PI=$((acos(-1.0)))

# http://maps.google.com.au/?ll=-32,116&z=16 Hartfield Park
# http://maps.google.com.au/?ll=-32,116&z=9 Metro
# http://maps.google.com.au/?ll=-30,120&z=6 SW
# http://maps.google.com.au/?ll=-25,135&z=4 AU

# ll resolution of 12' or 0.2°

# S = -Ve, N = +Ve
lat=$((-32.0 * $PI / 180.0))
coslat=$((cos($lat)))
sinlat=$((sin($lat)))

# W = -Ve, E = +Ve
long=$((116.0 * $PI / 180.0))

# TZ resolution of 10 minutes for SA and NT who are very special

# WA No daylight savings - ever
tz=8.0
# Very special people in CED and EUC
#tz=$((25.0/3.0))
# If there were normal people in SA and NT
#tz=9.0
# Special people in SA and NT
#tz=9.5
# Normal in QLD
#tz=10.0
# DS in SA
#tz=10.5
# DS NSW, VIC and TAS
#tz=11.0

tza=$(($tz * $PI / 12.0))

########################################################################
cat <<%%
Content-type: text/plain; charset=ISO-8859-1

In Perth according to the degree-day method (which neglects humidity):

 1. Cooling season is Nov-Apr 5 K·d at 24 °C where shading desired.

 2. Heating season is May-Oct 5 K·d at 12 °C where solar gain desired.

See <http://www.bom.gov.au/jsp/ncc/climate_averages/degree-days/>.
(It actually has the dimensions of temperature K·d/a or K·d/month
so it can be converted to cost and is the basis of AS 2627.1.)

The original Carrier method looked at reducing humidity. With such
latent cooling air-tightness is far more important than thermal
insulation, however after reducing humidity modern air conditioning then
goes on to provide sensible cooling (that is sensed with a dry bulb
thermostat). As we create a temperature difference heat will then flow
from hot to cold. This heat flow is only limited by thermal insulation.

When heating or cooling the inside temperature is held constant. As a
result of this energy consumption is proportional to the temperature
difference between inside and outside. It is the R value alone that
limits this rate of energy consumption. I.e. the degree-day method is
thus accurate. The (AS 2627.1) degree-day method only considers R value
it does not consider air infiltration, latent heat, heat capacities or
"thermal mass".

In passive design the inside temperature must vary. In this case latent
heat, the air tightness required to retain it, insulation levels and
even heat capacity are significant.

Some very simple rules on using windows need to be followed: --

 1. Cool in summer by opening windows at night.

 2. Heat in winter by having correctly shaded north facing windows.

In Perth metro area (-32,116) the sun is always to the north throughout
the year so south facing windows do not require any shade.

Windows facing to east or west will always get sun and must be avoided
or kept as small as possible. Such windows could benefit by having low
solar gain coatings and/or shading to reduce undesirable solar gain.
These windows may be required to provide day lighting or ventilation.

Windows facing to the north require shading as below to avoid summer
solar gain, but enable winter solar gain.

I.e. during the day: --

	-90° > azimuth < 90°

Solstices are close to days 173 and day 356.

	Sun direction at solar noon

	|        	|
	| 20°/70°	| 45°/45°
	|________	|________
	  Summer 	  Winter

Day 295,356,51 112,173,234 (see solar calcs below)
Mon Nov-Apr    May-Oct
HSA 70°        45°
E/S 0.36       1.00 (E/S should be between these values)

    E    Sun
\\_ |--|   /
~~\\_  |  /
=== \\_| /
 o H__\\/____
   |  /   |
 a | / h  S                     North --->
   |/_____|_ Top of sill
 A H ShadowAngle
   H
####--------

	          o
	tan(A) = ---, where A = 90° - ShadowAngle
	          a

Shading calculations with orientation error from true north to ±10°: --

%%
########################################################################

printf "Lat: % 5.1f° Long: % 5.1f° TZ: % 5.1fh standard time\n" \
$(($lat * 180.0 / $PI)) $(($long * 180.0 / $PI)) \
$(($tza * 12 / $PI))

# 1 Jan = 0
for (( day = 51.0 ; day <= 366.0 ; day += 61.0 ))
do

# Sun declination N = +Ve, S = -Ve approx
dec=$(((-23.44 * $PI / 180.0) * cos((2.0 * $PI / 365.0) * ($day + 10.0))))
sindec=$((sin($dec)))
cosdec=$((cos($dec)))

((isnan($dec))) && echo "Declination: $dec"

for (( orid = -10.0 ; orid <= 10.0; orid += 5.0 ))
do
ori=$(($orid * $PI / 180.0))

echo "\
D.ofYear Declination SolarTime  StandardT. ElevationA. AzimuthAngle ShadowAngle\
 Eave/Sill Orientation"

# sidereal time - solar time not TZ so long does not matter
#for hour in `seq 6 18`
for (( hour = 6.0 ; hour <= 18.0 ; hour += 2.0 ))
do

houra=$(((12 - $hour) * $PI / 12))

##houra=$(($houra + $tza - $long))

sinhoura=$((sin($houra)))

tzhoura=$(($houra + $long - $tza))

coshourdec=$((cos($houra) * $cosdec))

coshourdeclat=$(($coshourdec * $coslat))
sindeclat=$(($sindec * $sinlat))

# Solar elevation
sinelev=$(($coshourdeclat + $sindeclat))
elev=$((asin($sinelev)))

((isnan($elev))) && echo "Elevation: $elev"

(($elev < 0.0)) && continue

# Restructure this?
coselev=$((cos($elev)))
((isnan($coselev))) && echo "cos(elev): $coselev"
coselev2=$((sqrt(1 - ($sinelev * $sinelev))))
((isnan($coselev2))) && echo "cos(elev): $coselev2"
((abs($coselev - $coselev2) < 0.00001)) || \
{ print "coselev: $coselev = $coselev2 ?"; }

# Azimuth goes CW E of N +Ve W of N -Ve
# Use -Ve numbers for W of N to avoid discontinuity
# We only need one
nsinazim=$(($sinhoura * $cosdec / $coselev))
azim=$((asin($nsinazim)))

# Need 2 for N of the equator
#coshourdecsinlat=$(($coshourdec * $sinlat))
#sindeccoslat=$(($sindec * $coslat))
#cosazim=$((($coshourdecsinlat - $sindeccoslat) / $coselev))
#azim=$((acos($cosazim)))

((isnan($azim))) && echo "Azimuth: $azim"

tzhour=$((12 - ($tzhoura * 12 / $PI)))

# Eve to sill ratio
hsa=$(($azim - $ori))
tanvsa=$((tan($elev) / cos($hsa)))
vsa=$((atan($tanvsa)))
atanetos=$((($PI / 2) - $vsa))
etos=$((tan($atanetos)))

printf "Day:% 4.0f Dec:% 6.1f° Sol:% 5.1fh Std:% 5.1fh \
Elev:% 5.1f° Azim:% 6.1f° Shad:% 5.1f° E/S:% 5.2f Ori:% 5.1f°\n" \
$day $(($dec * 180.0 / $PI)) $hour $tzhour \
$(($elev * 180.0 / $PI)) $(($azim * 180.0 / $PI)) \
$(($vsa * 180.0 / $PI)) \
$etos $(($ori * 180.0 / $PI))

done
done
done