summaryrefslogtreecommitdiffstats
path: root/daylight.scm
diff options
context:
space:
mode:
Diffstat (limited to 'daylight.scm')
-rw-r--r--daylight.scm356
1 files changed, 356 insertions, 0 deletions
diff --git a/daylight.scm b/daylight.scm
new file mode 100644
index 0000000..6c989b2
--- /dev/null
+++ b/daylight.scm
@@ -0,0 +1,356 @@
+;;; "daylight.scm" Model of sun and sky colors.
+; Copyright 2001 Aubrey Jaffer
+;
+;Permission to copy this software, to modify it, to redistribute it,
+;to distribute modified versions, and to use it for any purpose is
+;granted, subject to the following restrictions and understandings.
+;
+;1. Any copy made of this software must include this copyright notice
+;in full.
+;
+;2. I have made no warranty or representation that the operation of
+;this software will be error-free, and I am under no obligation to
+;provide any services, by way of maintenance, update, or otherwise.
+;
+;3. In conjunction with products arising from the use of this
+;material, there shall be no use of my name in any advertising,
+;promotional, or sales literature without prior written consent in
+;each case.
+
+(require 'color-space)
+
+(define pi (* 4 (atan 1)))
+(define pi/180 (/ pi 180))
+
+;;@code{(require 'daylight)}
+;;@ftindex daylight
+;;@ftindex sunlight
+;;@ftindex sun
+;;@ftindex sky
+;;
+;;@noindent
+;;This package calculates the colors of sky as detailed in:@*
+;;@uref{http://www.cs.utah.edu/vissim/papers/sunsky/sunsky.pdf}@*
+;;@cite{A Practical Analytic Model for Daylight}@*
+;;A. J. Preetham, Peter Shirley, Brian Smits
+
+;;@body
+;;
+;;Returns the solar-time in hours given the integer @1 in the range 1 to
+;;366, and the local time in hours.
+;;
+;;To be meticulous, subtract 4 minutes for each degree of longitude west
+;;of the standard meridian of your time zone.
+(define (solar-hour julian-day hour)
+ (+ hour
+ (* 0.170 (sin (* 4 pi (- julian-day 80) 1/373)))
+ (* -0.129 (sin (* 2 pi (- julian-day 8) 1/355)))))
+
+;;@body
+(define (solar-declination julian-day)
+ (/ (* 0.4093 (sin (* 2 pi (- julian-day 81) 1/368))) pi/180))
+
+;;@body Returns a list of @var{theta_s}, the solar angle from the
+;;zenith, and @var{phi_s}, the solar azimuth. 0 <= @var{theta_s}
+;;measured in degrees. @var{phi_s} is measured in degrees from due
+;;south; west of south being positive.
+(define (solar-polar declination latitude solar-hour)
+ (define l (* pi/180 latitude))
+ (define d (* pi/180 declination))
+ (define pi*t/12 (* pi solar-hour 1/12))
+ (map (lambda (x) (/ x pi/180))
+ (list (- (/ pi 2) (asin (- (* (sin l) (sin d))
+ (* (cos l) (cos d) (cos pi*t/12)))))
+ (atan (* -1 (cos d) (sin pi*t/12))
+ (- (* (cos l) (sin d))
+ (* (sin l) (cos d) (cos pi*t/12)))))))
+
+;;@noindent
+;;In the following procedures, the number 0 <= @var{theta_s} <= 90 is
+;;the solar angle from the zenith in degrees.
+
+;;(plot (lambda (t) (+ -.5 (/ 9 (expt 1.55 t)))) 0 6) ;tweaked
+
+;;@cindex turbidity
+;;@noindent
+;;Turbidity is a measure of the fraction of scattering due to haze as
+;;opposed to molecules. This is a convenient quantity because it can be
+;;estimated based on visibility of distant objects. This model fails
+;;for turbidity values less than 1.3.
+;;
+;;@example
+;;@group
+;; _______________________________________________________________
+;;512|-: |
+;; | * pure-air |
+;;256|-:** |
+;; | : ** exceptionally-clear |
+;;128|-: * |
+;; | : ** |
+;; 64|-: * |
+;; | : ** very-clear |
+;; 32|-: ** |
+;; | : ** |
+;; 16|-: *** clear |
+;; | : **** |
+;; 8|-: **** |
+;; | : **** light-haze |
+;; 4|-: **** |
+;; | : ****** |
+;; 2|-: ******** haze thin-|
+;; | : *********** fog |
+;; 1|-:----------------------------------------------------*******--|
+;; |_:____.____:____.____:____.____:____.____:____.____:____.____:_|
+;; 1 2 4 8 16 32 64
+;; Meterorological range (km) versus Turbidity
+;;@end group
+;;@end example
+
+(define sol-spec
+ '#(16559.0
+ 16233.7
+ 21127.5
+ 25888.2
+ 25829.1
+ 24232.3
+ 26760.5
+ 29658.3
+ 30545.4
+ 30057.5
+ 30663.7
+ 28830.4
+ 28712.1
+ 27825.0
+ 27100.6
+ 27233.6
+ 26361.3
+ 25503.8
+ 25060.2
+ 25311.6
+ 25355.9
+ 25134.2
+ 24631.5
+ 24173.2
+ 23685.3
+ 23212.1
+ 22827.7
+ 22339.8
+ 21970.2
+ 21526.7
+ 21097.9
+ 20728.3
+ 20240.4
+ 19870.8
+ 19427.2
+ 19072.4
+ 18628.9
+ 18259.2
+ 17960 ;guesses for the rest
+ 17730
+ 17570))
+
+(define k_o-spec
+ '#(0.003
+ 0.006
+ 0.009
+ 0.014
+ 0.021
+ 0.03
+ 0.04
+ 0.048
+ 0.063
+ 0.075
+ 0.085
+ 0.103
+ 0.12
+ 0.12
+ 0.115
+ 0.125
+ 0.12
+ 0.105
+ 0.09
+ 0.079
+ 0.067
+ 0.057
+ 0.048
+ 0.036
+ 0.028
+ 0.023
+ 0.018
+ 0.014
+ 0.011
+ 0.01
+ 0.009
+ 0.007
+ 0.004
+ 0))
+
+;;@body Returns a vector of 41 values, the spectrum of sunlight from
+;;380.nm to 790.nm for a given @1 and @2.
+(define (sunlight-spectrum turbidity theta_s)
+ (define (solCurve wl) (vector-ref sol-spec (quotient (- wl 380) 10)))
+ (define (k_oCurve wl) (if (>= wl 450)
+ (vector-ref k_o-spec (quotient (- wl 450) 10))
+ 0))
+ (define (k_gCurve wl) (case wl
+ ((760) 3.0)
+ ((770) 0.21)
+ (else 0)))
+ (define (k_waCurve wl) (case wl
+ ((690) 0.016)
+ ((700) 0.024)
+ ((710) 0.0125)
+ ((720) 1)
+ ((730) 0.87)
+ ((740) 0.061)
+ ((750) 0.001)
+ ((760) 1.e-05)
+ ((770) 1.e-05)
+ ((780) 0.0006)
+ (else 0)))
+
+ (define data (make-vector (+ 1 (quotient (- 780 380) 10)) 0.0))
+ ;;alpha - ratio of small to large particle sizes. (0:4,usually 1.3)
+ (define alpha 1.3)
+ ;;beta - amount of aerosols present
+ (define beta (- (* 0.04608365822050 turbidity) 0.04586025928522))
+ ;;lOzone - amount of ozone in cm(NTP)
+ (define lOzone .35)
+ ;;w - precipitable water vapor in centimeters (standard = 2)
+ (define w 2.0)
+ ;;m - Relative Optical Mass
+ (define m (/ (+ (cos (* pi/180 theta_s))
+ (* 0.15 (expt (- 93.885 theta_s) -1.253)))))
+ (and
+ (not (negative? (- 93.885 theta_s)))
+ ;; Compute specturm of sunlight
+ (do ((wl 780 (+ -5 wl)))
+ ((< wl 380) data)
+ (let* (;;Rayleigh Scattering
+ ;; paper and program disagree!! Looks like font-size typo in paper.
+ ;;(tauR (exp (* -0.008735 (expt (/ wl 1000) (* -4.08 m))))) ;sunsky.pdf
+ (tauR (exp (* -0.008735 m (expt (/ wl 1000) -4.08)))) ;RiSunConstants.C
+ ;;Aerosal (water + dust) attenuation
+ ;; paper and program disagree!! Looks like font-size typo in paper.
+ ;;(tauA (exp (* -1 beta (expt (/ wl 1000) (* -1 m alpha)))))
+ (tauA (exp (* -1 m beta (expt (/ wl 1000) (- alpha)))))
+ ;;Attenuation due to ozone absorption
+ (tauO (exp (* -1 m (k_oCurve wl) lOzone)))
+ ;;Attenuation due to mixed gases absorption
+ (tauG (exp (* -1.41 m (k_gCurve wl)
+ (expt (+ 1 (* 118.93 m (k_gCurve wl))) -0.45))))
+ ;;Attenuation due to water vapor absorbtion
+ (tauWA (exp (* -0.2385 m w (k_waCurve wl)
+ (expt (+ 1 (* 20.07 m w (k_waCurve wl))) -0.45)))))
+ (vector-set! data (quotient (- wl 380) 10)
+ (* (solCurve wl) tauR tauA tauO tauG tauWA))))))
+
+;;@body Returns (unnormalized) XYZ values for color of sunlight for a
+;;given @1 and @2.
+(define (sunlight-XYZ turbidity theta_s)
+ (define spectrum (sunlight-spectrum turbidity theta_s))
+ (and spectrum (spectrum->XYZ spectrum 380.e-9 780.e-9)))
+
+;;@body Given @1 and @2, @0 returns the CIEXYZ triple for color of
+;;sunlight scaled to be just inside the RGB709 gamut.
+(define (sunlight-CIEXYZ turbidity theta_s)
+ (define spectrum (sunlight-spectrum turbidity theta_s))
+ (and spectrum (spectrum->CIEXYZ spectrum 380.e-9 780.e-9)))
+
+;; Arguments and result in radians
+(define (angle-between theta phi theta_s phi_s)
+ (define cospsi (+ (* (sin theta) (sin theta_s) (cos (- phi phi_s)))
+ (* (cos theta) (cos theta_s))))
+ (cond ((> cospsi 1) 0)
+ ((< cospsi -1) pi)
+ (else (acos cospsi))))
+
+;;@body Returns the xyY (chromaticity and luminance) at the zenith. The
+;;Luminance has units kcd/m^2.
+(define (zenith-xyY turbidity theta_s)
+ (let* ((ths (* theta_s pi/180))
+ (thetas (do ((th 1 (* ths th))
+ (lst '() (cons th lst))
+ (cnt 3 (+ -1 cnt)))
+ ((negative? cnt) lst)))
+ (turbds (do ((tr 1 (* turbidity tr))
+ (lst '() (cons tr lst))
+ (cnt 2 (+ -1 cnt)))
+ ((negative? cnt) lst))))
+ (append (map (lambda (row) (apply + (map * row turbds)))
+ (map color:linear-transform
+ '(((+0.00165 -0.00374 +0.00208 +0 )
+ (-0.02902 +0.06377 -0.03202 +0.00394)
+ (+0.11693 -0.21196 +0.06052 +0.25885))
+ ((+0.00275 -0.00610 +0.00316 +0 )
+ (-0.04214 +0.08970 -0.04153 +0.00515)
+ (+0.15346 -0.26756 +0.06669 +0.26688)))
+ (list thetas thetas)))
+ (list (+ (* (tan (* (+ 4/9 (/ turbidity -120)) (+ pi (* -2 ths))))
+ (- (* 4.0453 turbidity) 4.9710))
+ (* -0.2155 turbidity)
+ 2.4192)))))
+
+;;@body @1 is a positive real number expressing the amount of light
+;;scattering. The real number @2 is the solar angle from the zenith in
+;;degrees.
+;;
+;;@0 returns a function of one angle @var{theta}, the angle from the
+;;zenith of the viewing direction (in degrees); and returning the xyY
+;;value for light coming from that elevation of the sky.
+(define (overcast-sky-color-xyY turbidity theta_s)
+ (define xyY_z (zenith-xyY turbidity theta_s))
+ (lambda (theta . phi)
+ (list (car xyY_z) (cadr xyY_z)
+ (* 1/3 (caddr xyY_z) (+ 1 (* 2 (cos (* pi/180 theta))))))))
+
+;;@body @1 is a positive real number expressing the amount of light
+;;scattering. The real number @2 is the solar angle from the zenith in
+;;degrees. The real number @3 is the solar angle from south.
+;;
+;;@0 returns a function of two angles, @var{theta} and @var{phi} which
+;;specify the angles from the zenith and south meridian of the viewing
+;;direction (in degrees); returning the xyY value for light coming from
+;;that direction of the sky.
+;;
+;;@code{sky-color-xyY} calls @code{overcast-sky-color-xyY} for
+;;@1 <= 20; otherwise the @0 function.
+(define (clear-sky-color-xyY turbidity theta_s phi_s)
+ (define xyY_z (zenith-xyY turbidity theta_s))
+ (define th_s (* pi/180 theta_s))
+ (define ph_s (* pi/180 phi_s))
+ (define (F~ A B C D E)
+ (lambda (th gm)
+ (* (+ 1 (* A (exp (/ B (cos th)))))
+ (+ 1 (* C (exp (* D gm))) (* E (expt (cos gm) 2))))))
+ (let* ((tb1 (list turbidity 1))
+ (Fs (map (lambda (mat) (apply F~ (color:linear-transform mat tb1)))
+ '((( 0.17872 -1.46303)
+ (-0.35540 +0.42749)
+ (-0.02266 +5.32505)
+ ( 0.12064 -2.57705)
+ (-0.06696 +0.37027))
+ ((-0.01925 -0.25922)
+ (-0.06651 +0.00081)
+ (-0.00041 +0.21247)
+ (-0.06409 -0.89887)
+ (-0.00325 +0.04517))
+ ((-0.01669 -0.26078)
+ (-0.09495 +0.00921)
+ (-0.00792 +0.21023)
+ (-0.04405 -1.65369)
+ (-0.01092 +0.05291)))))
+ (F_0s (map (lambda (F) (F 0 th_s)) Fs)))
+ (lambda (theta phi)
+ (let* ((th (* pi/180 theta))
+ (ph (* pi/180 phi))
+ (gm (angle-between th_s ph_s th ph)))
+ ;;(print th ph '=> gm)
+ (map (lambda (x F F_0) (* x (/ (F th gm) F_0)))
+ xyY_z
+ Fs
+ F_0s)))))
+(define (sky-color-xyY turbidity theta_s phi_s)
+ (if (> turbidity 20)
+ (overcast-sky-color-xyY turbidity theta_s)
+ (clear-sky-color-xyY turbidity theta_s phi_s)))