diff options
Diffstat (limited to 'daylight.scm')
| -rw-r--r-- | daylight.scm | 356 |
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))) |