| 1 | 
greg | 
2.1 | 
#!/usr/bin/perl -w | 
| 2 | 
greg | 
2.5 | 
# RCSid $Id: genskyvec.pl,v 2.4 2010/01/23 22:30:29 greg Exp $ | 
| 3 | 
greg | 
2.1 | 
# | 
| 4 | 
  | 
  | 
# Generate Reinhart vector for a given sky description | 
| 5 | 
  | 
  | 
# | 
| 6 | 
  | 
  | 
#       G. Ward | 
| 7 | 
  | 
  | 
# | 
| 8 | 
  | 
  | 
use strict; | 
| 9 | 
  | 
  | 
my @skycolor = (0.960, 1.004, 1.118); | 
| 10 | 
greg | 
2.3 | 
my $mf = 4; | 
| 11 | 
greg | 
2.1 | 
while ($#ARGV >= 0) { | 
| 12 | 
  | 
  | 
        if ("$ARGV[0]" eq "-c") { | 
| 13 | 
  | 
  | 
                @skycolor = @ARGV[1..3]; | 
| 14 | 
  | 
  | 
                shift @ARGV; shift @ARGV; shift @ARGV; | 
| 15 | 
greg | 
2.3 | 
        } elsif ("$ARGV[0]" eq "-m") { | 
| 16 | 
  | 
  | 
                $mf = $ARGV[1]; | 
| 17 | 
  | 
  | 
                shift @ARGV; | 
| 18 | 
greg | 
2.1 | 
        } | 
| 19 | 
  | 
  | 
        shift @ARGV; | 
| 20 | 
  | 
  | 
} | 
| 21 | 
  | 
  | 
# Load sky description into line array, separating sun if one | 
| 22 | 
  | 
  | 
my @skydesc; | 
| 23 | 
  | 
  | 
my $lightline; | 
| 24 | 
  | 
  | 
my @sunval; | 
| 25 | 
  | 
  | 
my $sunline; | 
| 26 | 
  | 
  | 
my $skyOK = 0; | 
| 27 | 
  | 
  | 
my $srcmod;     # putting this inside loop breaks code(?!) | 
| 28 | 
  | 
  | 
while (<>) { | 
| 29 | 
  | 
  | 
        push @skydesc, $_; | 
| 30 | 
  | 
  | 
        if (/^\w+\s+light\s+/) { | 
| 31 | 
  | 
  | 
                s/\s*$//; s/^.*\s//; | 
| 32 | 
  | 
  | 
                $srcmod = $_; | 
| 33 | 
  | 
  | 
                $lightline = $#skydesc; | 
| 34 | 
  | 
  | 
        } elsif (defined($srcmod) && /^($srcmod)\s+source\s/) { | 
| 35 | 
greg | 
2.5 | 
                @sunval = split(' ', $skydesc[$lightline + 3]); | 
| 36 | 
greg | 
2.1 | 
                shift @sunval; | 
| 37 | 
  | 
  | 
                $sunline = $#skydesc; | 
| 38 | 
  | 
  | 
        } elsif (/\sskyfunc\s*$/) { | 
| 39 | 
  | 
  | 
                $skyOK = 1; | 
| 40 | 
  | 
  | 
        } | 
| 41 | 
  | 
  | 
} | 
| 42 | 
  | 
  | 
die "Bad sky description!\n" if (! $skyOK); | 
| 43 | 
  | 
  | 
# Strip out the solar source if present | 
| 44 | 
  | 
  | 
my @sundir; | 
| 45 | 
  | 
  | 
if (defined $sunline) { | 
| 46 | 
greg | 
2.5 | 
        @sundir = split(' ', $skydesc[$sunline + 3]); | 
| 47 | 
greg | 
2.1 | 
        shift @sundir; | 
| 48 | 
  | 
  | 
        undef @sundir if ($sundir[2] <= 0); | 
| 49 | 
  | 
  | 
        splice(@skydesc, $sunline, 5); | 
| 50 | 
  | 
  | 
} | 
| 51 | 
  | 
  | 
# Reinhart sky sample generator | 
| 52 | 
  | 
  | 
my $rhcal = ' | 
| 53 | 
  | 
  | 
DEGREE : PI/180; | 
| 54 | 
  | 
  | 
x1 = .5; x2 = .5; | 
| 55 | 
  | 
  | 
alpha : 90/(MF*7 + .5); | 
| 56 | 
  | 
  | 
tnaz(r) : select(r, 30, 30, 24, 24, 18, 12, 6); | 
| 57 | 
  | 
  | 
rnaz(r) : if(r-(7*MF-.5), 1, MF*tnaz(floor((r+.5)/MF) + 1)); | 
| 58 | 
  | 
  | 
raccum(r) : if(r-.5, rnaz(r-1) + raccum(r-1), 0); | 
| 59 | 
  | 
  | 
RowMax : 7*MF + 1; | 
| 60 | 
  | 
  | 
Rmax : raccum(RowMax); | 
| 61 | 
  | 
  | 
Rfindrow(r, rem) : if(rem-rnaz(r)-.5, Rfindrow(r+1, rem-rnaz(r)), r); | 
| 62 | 
greg | 
2.2 | 
Rrow = if(Rbin-(Rmax-.5), RowMax-1, Rfindrow(0, Rbin)); | 
| 63 | 
greg | 
2.1 | 
Rcol = Rbin - raccum(Rrow) - 1; | 
| 64 | 
  | 
  | 
Razi_width = 2*PI / rnaz(Rrow); | 
| 65 | 
  | 
  | 
RAH : alpha*DEGREE; | 
| 66 | 
  | 
  | 
Razi = if(Rbin-.5, (Rcol + x2 - .5)*Razi_width, 2*PI*x2); | 
| 67 | 
greg | 
2.2 | 
Ralt = if(Rbin-.5, (Rrow + x1)*RAH, asin(-x1)); | 
| 68 | 
greg | 
2.1 | 
Romega = if(.5-Rbin, 2*PI, if(Rmax-.5-Rbin,  | 
| 69 | 
  | 
  | 
        Razi_width*(sin(RAH*(Rrow+1)) - sin(RAH*Rrow)), | 
| 70 | 
  | 
  | 
        2*PI*(1 - cos(RAH/2)) ) ); | 
| 71 | 
  | 
  | 
cos_ralt = cos(Ralt); | 
| 72 | 
  | 
  | 
Dx = sin(Razi)*cos_ralt; | 
| 73 | 
  | 
  | 
Dy = cos(Razi)*cos_ralt; | 
| 74 | 
  | 
  | 
Dz = sin(Ralt); | 
| 75 | 
  | 
  | 
'; | 
| 76 | 
greg | 
2.3 | 
my $nbins = `rcalc -n -e MF:$mf -e \'$rhcal\' -e \'\$1=Rmax+1\'`; | 
| 77 | 
  | 
  | 
chomp $nbins; | 
| 78 | 
greg | 
2.1 | 
# Create octree for rtrace | 
| 79 | 
  | 
  | 
my $octree = "/tmp/gtv$$.oct"; | 
| 80 | 
  | 
  | 
open OCONV, "| oconv - > $octree"; | 
| 81 | 
  | 
  | 
print OCONV @skydesc; | 
| 82 | 
  | 
  | 
print OCONV "skyfunc glow skyglow 0 0 4 @skycolor 0\n"; | 
| 83 | 
  | 
  | 
print OCONV "skyglow source sky 0 0 4 0 0 1 360\n"; | 
| 84 | 
  | 
  | 
close OCONV; | 
| 85 | 
  | 
  | 
# Run rtrace and average output for every 16 samples | 
| 86 | 
greg | 
2.3 | 
my $tregcommand = "cnt $nbins 16 | rcalc -of -e MF:$mf -e '$rhcal' " . | 
| 87 | 
greg | 
2.1 | 
        q{-e 'Rbin=$1;x1=rand(recno*.37-5.3);x2=rand(recno*-1.47+.86)' } . | 
| 88 | 
  | 
  | 
        q{-e '$1=0;$2=0;$3=0;$4=Dx;$5=Dy;$6=Dz' } . | 
| 89 | 
  | 
  | 
        "| rtrace -h -ff -ab 0 -w $octree | total -if3 -16 -m"; | 
| 90 | 
  | 
  | 
my @tregval = `$tregcommand`; | 
| 91 | 
  | 
  | 
unlink $octree; | 
| 92 | 
  | 
  | 
# Find closest 3 patches to sun and divvy up direct solar contribution | 
| 93 | 
  | 
  | 
my @bestdir; | 
| 94 | 
  | 
  | 
if (@sundir) { | 
| 95 | 
  | 
  | 
        my $somega = ($sundir[3]/360)**2 * 3.141592654**3; | 
| 96 | 
greg | 
2.2 | 
        my $cmd = "cnt " . ($nbins-1) . | 
| 97 | 
greg | 
2.3 | 
                " | rcalc -e MF:$mf -e '$rhcal' -e Rbin=recno " . | 
| 98 | 
greg | 
2.1 | 
                "-e 'dot=Dx*$sundir[0] + Dy*$sundir[1] + Dz*$sundir[2]' " . | 
| 99 | 
  | 
  | 
                "-e 'cond=dot-.866' " . | 
| 100 | 
  | 
  | 
                q{-e '$1=if(1-dot,acos(dot),0);$2=Romega;$3=recno' }; | 
| 101 | 
  | 
  | 
        @bestdir = `$cmd | sort -n | head -3`; | 
| 102 | 
  | 
  | 
        my (@ang, @dom, @ndx); | 
| 103 | 
  | 
  | 
        my $wtot = 0; | 
| 104 | 
  | 
  | 
        for my $i (0..2) { | 
| 105 | 
greg | 
2.5 | 
                ($ang[$i], $dom[$i], $ndx[$i]) = split(' ', $bestdir[$i]); | 
| 106 | 
greg | 
2.1 | 
                $wtot += 1./($ang[$i]+.02); | 
| 107 | 
  | 
  | 
        } | 
| 108 | 
  | 
  | 
        for my $i (0..2) { | 
| 109 | 
  | 
  | 
                my $wt = 1./($ang[$i]+.02)/$wtot * $somega / $dom[$i]; | 
| 110 | 
greg | 
2.5 | 
                my @scolor = split(' ', $tregval[$ndx[$i]]); | 
| 111 | 
greg | 
2.1 | 
                for my $j (0..2) { $scolor[$j] += $wt * $sunval[$j]; } | 
| 112 | 
greg | 
2.4 | 
                $tregval[$ndx[$i]] = "$scolor[0]\t$scolor[1]\t$scolor[2]\n"; | 
| 113 | 
greg | 
2.1 | 
        } | 
| 114 | 
  | 
  | 
} | 
| 115 | 
  | 
  | 
# Output our final vector | 
| 116 | 
  | 
  | 
print @tregval; |