Sub Command1_Click ()
' PROGRAM VB_SKY_R  Sky Color with Clouds, ozone, volcano + aerosols, Eclipse
' VERSION BEGUN 21 MAY 2004 COMPLETED 07 SEPT 2004 BY STANLEY DAVID GEDZELMAN
' REFRACTION ADDED 15-23 OCT, CLOUDS CORRECTED 30 OCT 2004   HOORAY!!!
' PARAMETERIZATIONS FOR REFRACTION IN EXCEL FILES SBEAMTAU, OBS_DH AND PULKV
' CLOUD SHADOWS ADDED 14 DEC 2004 SKY VALID BELOW PHICLD
'*********** FUNCTION TANH  put in GENERAL
'Function tanh(x as single) as single
'tanh = (Exp(x) - Exp(-x)) / (Exp(x) + Exp(-x))
'End function
' **************************

Const TK0 = 273
Dim Kray As Single, p0 As Single, rd As Single
Dim p As Single, pold As Single
Dim T0 As Single, Trat As Single
Dim KAer As Single, beta As Single, Kvul As Single
Dim n0 As Single, nrf As Single, rfi As Single
Dim rho0 As Single, g0 As Single, Re As Single, rec As Single
Dim Re0 As Single, Retop As Single
Dim Hscl As Single, Hpscl As Single, HAer As Single, HVul As Single
Dim h0 As Single, hV0 As Single
Dim tau As Single, dtauin As Single, dtauout As Single
Dim TauA As Single, taur As Single
Dim tauR0 As Single, TauA0 As Single
Dim tauVul As Single, dtauV As Single, dtauZ As Single
Dim tauVSun As Single, TauZsun As Single
Dim dtauVc As Single, dtauVcin As Single
Dim tauRrat As Single, tauArat As Single
Dim zReff As Single, zAeff As Single
Dim kLV As Single, kLA As Single, kls As Single
Dim pi As Single, drc As Single
Dim xp As Single, xp0 As Single
Dim dxp As Single, dsp As Single, dx2 As Single
Dim hp As Single, dhp As Single
Dim zp As Single, zp0 As Single
Dim xp2 As Single, hp2 As Single, cc As Single
Dim xf As Single, zf As Single, hf As Single
Dim xs As Single, xa As Single
Dim xt As Single, yt As Single, zt As Single
Dim xsun As Single, ysun As Single, zsun As Single
Dim xsunchrm As Single, ysunchrm As Single, Fsun As Single
Dim xrel As Single, yrel As Single, xpen As Single
Dim htop As Single, xtop As Single, dtht As Single
Dim sgRAY As Single, sgAER As Single, Fmax As Single

Dim r As Single, g As Single, B As Single
Dim rp As Single, xcolor As Single, ycolor As Single
Dim i As Integer, il As Integer, ith As Integer, ithmax As Integer
Dim ix As Integer, ix2 As Integer, ixmax As Integer, j As Integer
Dim ixpen As Integer, eclfac As Single
Dim eye As Single
Dim Fnorm As Single
ReDim lmd(61) As Single
ReDim lmdR(61) As Single, lmdA(61) As Single, lmdV(61) As Single
ReDim F0(61) As Single, F(61) As Single
ReDim irrad(100) As Single
ReDim xchrm(100) As Single, ychrm(100) As Single
Dim hZmx As Single, htrop As Single, HZ As Single                              '**
Dim nZ As Single, nZ2 As Single, nZmx As Single, nZtrop As Single
Dim sigmZ As Single, kZ As Single, lmdZ As Single, Dobson As Single
ReDim sgmz(100) As Single
ReDim hfp(101, 43) As Single, dxf(43) As Single
Dim mf As Single, mfx As Single
Const xfrm = 1350
' physical chromaticity coordinates
Const xbm = .1703, ybm = .0058
Const xgm = .1547, ygm = .8059
Const xrm = .726, yrm = .274
Const xneut = (xrm + xgm + xbm) / 3#
Const yneut = (yrm + ygm + ybm) / 3#
Dim aa As String' To read first line of chromaticity files
Dim thcie As Single
Dim detx As Single
Dim detR As Single, detG As Single, detB As Single
ReDim xv(61) As Single, yv(61) As Single, thv(61) As Single
ReDim xbar(61) As Single, ybar(61) As Single, zbar(61) As Single
Dim mp As Single, bp As Single
Dim pur As Single, purpct As Single
Dim xpur As Single, ypur As Single
Dim huemax As Single
' Sky beam variables
Dim iphi As Integer, iphimax As Integer
Dim tht0 As Single, tht00 As Single, dth As Single
Dim thtp As Single, thtcor As Single, SKYFRAC As Single
Dim zss As Single, thtss As Single
Dim dphi As Single, delphi As Single, phif As Single
Dim psi As Single, aAer As Single, aVul As Single
Dim bAer As Single, bVul As Single
Dim cA As Single, cV As Single
Dim x0 As Single, z0 As Single, x_ab As Single
Dim skyR As Single, skyA As Single, skyV As Single, skyZ As Single
Dim skyRc As Single, skyAc As Single, skyVc As Single
Dim skyRcin As Single, skyAcin As Single
Dim Fsct As Single
Dim czdraw As Single, rsun As Integer, rcld As Integer
ReDim phi(80) As Single, Fskychrm(61, 200) As Single
ReDim Fskytot(200) As Single
ReDim red(200), grn(200), blu(200)

' Cloud Variables And Constants
Dim hcld As Single, phicld As Single, zcloud As Single
Dim ixcld As Integer, iphicld As Integer
Dim ncld As Integer, ctr As Integer, ntrop As Integer
Dim nu As Integer, nd As Integer
Dim dxcld As Single, phic As Single, phicb As Single
Dim taucld As Single, taub As Single
Dim acld As Single, Pcld As Single, qq As Single
Dim albC As Single, trnC As Single, pwr As Single
Dim cldil As Single, edgil As Single
Dim albE As Single, trnE As Single
Dim dalb As Single, trnpC As Single, trnpE As Single
Dim xcld As Single, ycld As Single, zcld As Single
Dim xcldchrm As Single, ycldchrm As Single
Dim Fcld As Single, FMax2 As Single
Dim xedg As Single, yedg As Single, zedg As Single
Dim xedgchrm As Single, yedgchrm As Single, Fedg As Single
Dim thtcld As Single, IrrCld As Single, IrrEdg As Single, IrrSun As Single
ReDim Fcldchrm(61) As Single, Fedgchrm(61) As Single
ReDim alb(2) As Single

Randomize Timer
pi = 4 * Atn(1): drc = pi / 180#
rcld = 300
rsun = 100
Dim CHRM As Integer, CLOUD As Integer, REFRACT As Integer
Dim SHADOW As Integer, ECLIPSE As Integer, OZONE As Integer
Dim MULT As Integer
If option1.Value = True Then CHRM = 0 Else CHRM = 1
If CHECK2.Value = 0 Then CLOUD = 0 Else CLOUD = 1
If CHECK3.Value = 0 Then ECLIPSE = 0 Else ECLIPSE = 1
If CHECK4.Value = 0 Then MULT = 0 Else MULT = 1
If CHECK5.Value = 0 Then OZONE = 0 Else OZONE = 1
If CHECK6.Value = 0 Then REFRACT = 0 Else REFRACT = 1
If CHECK7.Value = 0 Then SHADOW = 0 Else SHADOW = 1
'*********************************************************
'  CHOOSE BASIC VARIABLES  +   CLOUDS                  '**
ixmax = 100' number of cells                           '**
iphimax = 42' number of viewing angles                 '**
SKYFRAC = hscroll12.Value / 20' Fraction of sky        '**
ithmax = 8'  number of solar angles                    '**
tht00 = hscroll1.Value + .01' START SOLAR ELEVATION ANGLE'**
dth = hscroll2.Value / 10'4                            '**
h0 = hscroll3.Value * 100' observer's height msl (m)   '**
If CHRM = 1 Then eye = 0 Else eye = .1                 '**
phicld = hscroll5.Value * drc ' CLOUD ELEVATION ANGLE  '**
hcld = hscroll4.Value * 100 + h0' CLOUD HT ABOVE VIEWER'**
taucld = hscroll6.Value ' CLOUD OPTICAL THICKNESS      '**
T0 = hscroll13.Value'             TEMPERATURE                      '**
p0 = 101325: rd = 287                                  '**
xpen = 80000# 'Distance to penumbra -center of eclipse '**
HAer = hscroll7.Value * 100' 3000 AEROSOL SCALE HEIGHT '**
'*********************************************************
mf = 7800: mfx = .8 / (1 - .2 * CHRM)
'drawwidth = 4
scaleheight = -mf / (1 + .333 * CHRM)'+ 600'+1200
scaletop = mf + 250 * (10 - 9 * CHRM)
scaleleft = -750
windowstate = 2
Line (-1000, 0)-(15000, 12000), RGB(255, 255, 255), BF

option1.Visible = False
option2.Visible = False
CHECK2.Visible = False
CHECK3.Visible = False
CHECK4.Visible = False
CHECK5.Visible = False
CHECK6.Visible = False
CHECK7.Visible = False
hscroll1.Visible = False
hscroll2.Visible = False
hscroll3.Visible = False
hscroll4.Visible = False
hscroll5.Visible = False
hscroll6.Visible = False
hscroll7.Visible = False
hscroll8.Visible = False
hscroll9.Visible = False
hscroll10.Visible = False
hscroll11.Visible = False
hscroll12.Visible = False
hscroll13.Visible = False
command1.Visible = False

Trat = TK0 / T0
rho0 = p0 / (rd * T0)
g0 = 9.81
Re = 6370000!: Rem effective radius should include refraction
' 80 km max height for beam - must be larger if phi near pi and tht0 near -8 deg
htop = 80000
Re0 = Re + h0
Retop = Re + htop
rec = 100
Hscl = rd * T0 / g0
Hpscl = 2 * Hscl * p0
Kray = 2 / 3 * (.00000083839)
n0 = .000293 * TK0 / T0
If REFRACT = 1 Then thtss = -.52 * Trat ^ 1.7 * drc Else thtss = 0' sunset point and refraction angle
acld = 3 ' angular focus factor for cloud droplets - min at pi
Pcld = 1 - Exp(-pi * acld)
ncld = 30000
phi(1) = -Sqr(1.5 * h0 / Re) + .001' first viewing angle above horizon
delphi = pi - 2 * phi(1) + .001
For iphi = 2 To iphimax
dphi = (.008 + SKYFRAC * .12 * Abs(Cos((phi(iphi - 1) + dphi / 2 - pi / 2) * pi / delphi)) ^ .6) * delphi / pi
phi(iphi) = phi(iphi - 1) + dphi
If Abs(phi(iphi) - phicld) < dphi / 1.8 Then iphicld = iphi
Next iphi
czdraw = 1.25 * mf / phi(iphimax)
' ********  CALCULATING CLOUD AND VIEWER ALBEDO *******
If CLOUD = 1 Then
For ix = 1 To 2
phicb = 1.56 * (ix - 1)' prepare to interpolate to tht0
For i = 1 To ncld
phic = phicb
ctr = 0
taub = 0
Do Until ctr > 0
qq = .00001 + .99998 * Rnd
taub = taub - Log(qq) * Cos(phic)
If taub > taucld Then
ctr = 1
nd = nd + 1
ElseIf taub < 0 Then
ctr = 2
nu = nu + 1
Else
qq = .00001 + .99998 * Rnd
phic = phic + Sgn(Rnd - .5) * Log(1 - qq * Pcld) / acld
End If
Loop
Next i
alb(ix) = nu / ncld
nd = 0: nu = 0
Next ix
pwr = 2 + Exp(tau - 1)
dalb = (alb(2) - alb(1)) / (1.56 ^ pwr)
trnpC = 1 - (alb(1) + dalb * Abs(1.56 - phicld) ^ pwr)
trnpE = 1 - .2 * (alb(1) + dalb * Abs(1.56 - phicld) ^ pwr)' Assumed fraction of Cloud
End If
lmd(1) = .4: xv(1) = .1733: yv(1) = .0048
lmd(2) = .405: xv(2) = .173: yv(2) = .00478
lmd(3) = .41: xv(3) = .1726: yv(3) = .0048
lmd(4) = .415: xv(4) = .1721: yv(4) = .0048
lmd(5) = .42: xv(5) = .1714: yv(5) = .0051
lmd(6) = .425: xv(6) = .1703: yv(6) = .0058
lmd(7) = .43: xv(7) = .1689: yv(7) = .0069
lmd(8) = .435: xv(8) = .1669: yv(8) = .0086
lmd(9) = .44: xv(9) = .1644: yv(9) = .0109
lmd(10) = .445: xv(10) = .1611: yv(10) = .0138
lmd(11) = .45: xv(11) = .1566: yv(11) = .0177
lmd(12) = .455: xv(12) = .151: yv(12) = .0227
lmd(13) = .46: xv(13) = .144: yv(13) = .0297
lmd(14) = .465: xv(14) = .1355: yv(14) = .0399
lmd(15) = .47: xv(15) = .1241: yv(15) = .0578
lmd(16) = .475: xv(16) = .1096: yv(16) = .0868
lmd(17) = .48: xv(17) = .0913: yv(17) = .1327
lmd(18) = .485: xv(18) = .0687: yv(18) = .2007
lmd(19) = .49: xv(19) = .0454: yv(19) = .295
lmd(20) = .495: xv(20) = .0235: yv(20) = .4127
lmd(21) = .5: xv(21) = .0082: yv(21) = .5384
lmd(22) = .505: xv(22) = .0039: yv(22) = .6548
lmd(23) = .51: xv(23) = .0139: yv(23) = .7502
lmd(24) = .515: xv(24) = .0389: yv(24) = .812
lmd(25) = .52: xv(25) = .0743: yv(25) = .8338
lmd(26) = .525: xv(26) = .1142: yv(26) = .8262
lmd(27) = .53: xv(27) = .1547: yv(27) = .8059
lmd(28) = .535: xv(28) = .1929: yv(28) = .7816
lmd(29) = .54: xv(29) = .2296: yv(29) = .7543
lmd(30) = .545: xv(30) = .2658: yv(30) = .7243
lmd(31) = .55: xv(31) = .3016: yv(31) = .6923
lmd(32) = .555: xv(32) = .3373: yv(32) = .6589
lmd(33) = .56: xv(33) = .3731: yv(33) = .6245
lmd(34) = .565: xv(34) = .4078: yv(34) = .5897
lmd(35) = .57: xv(35) = .4441: yv(35) = .5547
lmd(36) = .575: xv(36) = .4788: yv(36) = .5202
lmd(37) = .58: xv(37) = .5125: yv(37) = .4866
lmd(38) = .585: xv(38) = .5448: yv(38) = .4545
lmd(39) = .59: xv(39) = .5752: yv(39) = .4242
lmd(40) = .595: xv(40) = .6029: yv(40) = .3966
lmd(41) = .6: xv(41) = .627: yv(41) = .3726
lmd(42) = .605: xv(42) = .6482: yv(42) = .3515
lmd(43) = .61: xv(43) = .6658: yv(43) = .334
lmd(44) = .615: xv(44) = .6801: yv(44) = .3197
lmd(45) = .62: xv(45) = .6915: yv(45) = .3083
lmd(46) = .625: xv(46) = .7006: yv(46) = .2993
lmd(47) = .63: xv(47) = .7079: yv(47) = .292
lmd(48) = .635: xv(48) = .714: yv(48) = .2859
lmd(49) = .64: xv(49) = .719: yv(49) = .281
lmd(50) = .645: xv(50) = .723: yv(50) = .277
lmd(51) = .65: xv(51) = .726: yv(51) = .274
lmd(52) = .655: xv(52) = .7283: yv(52) = .2717
lmd(53) = .66: xv(53) = .73: yv(53) = .27
lmd(54) = .665: xv(54) = .7311: yv(54) = .2689
lmd(55) = .67: xv(55) = .732: yv(55) = .268
lmd(56) = .675: xv(56) = .7327: yv(56) = .2673
lmd(57) = .68: xv(57) = .7334: yv(57) = .2666
lmd(58) = .685: xv(58) = .734: yv(58) = .266
lmd(59) = .69: xv(59) = .7344: yv(59) = .2656
lmd(60) = .695: xv(60) = .7346: yv(60) = .2654
lmd(61) = .7: xv(61) = .7347: yv(61) = .2653

For il = 1 To 61
'Planck's Law for Sun at 5700 K
F0(il) = 2000# / (lmd(il) ^ 5# * (Exp(2.481 / lmd(il)) - 1#))
Next il
If CHRM = 1 Then
   mp = (yv(61) - yv(1)) / (xv(61) - xv(1))
   bp = yv(61) - mp * xv(61)
End If
lmd(1) = .4: xbar(1) = .0144: ybar(1) = .0004: zbar(1) = .0685
lmd(2) = .405: xbar(2) = .0233: ybar(2) = .00064: zbar(2) = .1105
lmd(3) = .41: xbar(3) = .0432: ybar(3) = .0012: zbar(3) = .2056
lmd(4) = .415: xbar(4) = .078: ybar(4) = .00218: zbar(4) = .373
lmd(5) = .42: xbar(5) = .1344: ybar(5) = .004: zbar(5) = .6459
lmd(6) = .425: xbar(6) = .2132: ybar(6) = .0073: zbar(6) = 1.0317
lmd(7) = .43: xbar(7) = .2839: ybar(7) = .0116: zbar(7) = 1.3856
lmd(8) = .435: xbar(8) = .3268: ybar(8) = .0168: zbar(8) = 1.6142
lmd(9) = .44: xbar(9) = .3469: ybar(9) = .023: zbar(9) = 1.7402
lmd(10) = .445: xbar(10) = .3483: ybar(10) = .0298: zbar(10) = 1.784
lmd(11) = .45: xbar(11) = .3326: ybar(11) = .038: zbar(11) = 1.7727
lmd(12) = .455: xbar(12) = .3193: ybar(12) = .048: zbar(12) = 1.7472
lmd(13) = .46: xbar(13) = .2909: ybar(13) = .06: zbar(13) = 1.6693
lmd(14) = .465: xbar(14) = .2509: ybar(14) = .0739: zbar(14) = 1.5268
lmd(15) = .47: xbar(15) = .1954: ybar(15) = .091: zbar(15) = 1.288
lmd(16) = .475: xbar(16) = .1422: ybar(16) = .1126: zbar(16) = 1.0427
lmd(17) = .48: xbar(17) = .0956: ybar(17) = .139: zbar(17) = .8128
lmd(18) = .485: xbar(18) = .058: ybar(18) = .1693: zbar(18) = .6163
lmd(19) = .49: xbar(19) = .032: ybar(19) = .208: zbar(19) = .4651
lmd(20) = .495: xbar(20) = .0147: ybar(20) = .2586: zbar(20) = .3532
lmd(21) = .5: xbar(21) = .0049: ybar(21) = .3232: zbar(21) = .272
lmd(22) = .505: xbar(22) = .0024: ybar(22) = .4073: zbar(22) = .2123
lmd(23) = .51: xbar(23) = .0093: ybar(23) = .503: zbar(23) = .1582
lmd(24) = .515: xbar(24) = .0291: ybar(24) = .6082: zbar(24) = .1117
lmd(25) = .52: xbar(25) = .0633: ybar(25) = .71: zbar(25) = .0783
lmd(26) = .525: xbar(26) = .1096: ybar(26) = .7932: zbar(26) = .0572
lmd(27) = .53: xbar(27) = .1655: ybar(27) = .8682: zbar(27) = .0421
lmd(28) = .535: xbar(28) = .2258: ybar(28) = .9194: zbar(28) = .0299
lmd(29) = .54: xbar(29) = .2904: ybar(29) = .954: zbar(29) = .0204
lmd(30) = .545: xbar(30) = .3597: ybar(30) = .9802: zbar(30) = .0134
lmd(31) = .55: xbar(31) = .4335: ybar(31) = .995: zbar(31) = .00877
lmd(32) = .555: xbar(32) = .512: ybar(32) = 1.0002: zbar(32) = .00577
lmd(33) = .56: xbar(33) = .5945: ybar(33) = .995: zbar(33) = .00386
lmd(34) = .565: xbar(34) = .6783: ybar(34) = .9786: zbar(34) = .00269
lmd(35) = .57: xbar(35) = .7622: ybar(35) = .952: zbar(35) = .00206
lmd(36) = .575: xbar(36) = .8425: ybar(36) = .9154: zbar(36) = .00172
lmd(37) = .58: xbar(37) = .9162: ybar(37) = .87: zbar(37) = .0066
lmd(38) = .585: xbar(38) = .9785: ybar(38) = .8163: zbar(38) = .00135
lmd(39) = .59: xbar(39) = 1.0266: ybar(39) = .757: zbar(39) = .00116
lmd(40) = .595: xbar(40) = 1.0566: ybar(40) = .6949: zbar(40) = .00096
lmd(41) = .6: xbar(41) = 1.062: ybar(41) = .6316: zbar(41) = .00076
lmd(42) = .605: xbar(42) = 1.0453: ybar(42) = .5668: zbar(42) = .00056
lmd(43) = .61: xbar(43) = 1.0028: ybar(43) = .503: zbar(43) = .00038
lmd(44) = .615: xbar(44) = .9387: ybar(44) = .4412: zbar(44) = .00028
lmd(45) = .62: xbar(45) = .8545: ybar(45) = .381: zbar(45) = .0002
lmd(46) = .625: xbar(46) = .7515: ybar(46) = .321: zbar(46) = .00014
lmd(47) = .63: xbar(47) = .6424: ybar(47) = .265: zbar(47) = .00009
lmd(48) = .635: xbar(48) = .5419: ybar(48) = .217: zbar(48) = .00005
lmd(49) = .64: xbar(49) = .4479: ybar(49) = .175: zbar(49) = .00003
lmd(50) = .645: xbar(50) = .3609: ybar(50) = .1382: zbar(50) = .00002
lmd(51) = .65: xbar(51) = .2835: ybar(51) = .107: zbar(51) = .00001
lmd(52) = .655: xbar(52) = .2186: ybar(52) = .0816: zbar(52) = 0
lmd(53) = .66: xbar(53) = .1649: ybar(53) = .061: zbar(53) = 0
lmd(54) = .665: xbar(54) = .1212: ybar(54) = .0446: zbar(54) = 0
lmd(55) = .67: xbar(55) = .0874: ybar(55) = .032: zbar(55) = 0
lmd(56) = .675: xbar(56) = .0637: ybar(56) = .0232: zbar(56) = 0
lmd(57) = .68: xbar(57) = .0468: ybar(57) = .017: zbar(57) = 0
lmd(58) = .685: xbar(58) = .0329: ybar(58) = .0119: zbar(58) = 0
lmd(59) = .69: xbar(59) = .0227: ybar(59) = .0082: zbar(59) = 0
lmd(60) = .695: xbar(60) = .0158: ybar(60) = .00572: zbar(60) = 0
lmd(61) = .7: xbar(61) = .0114: ybar(61) = .0041: zbar(61) = 0

detx = xrm * ygm + xgm * ybm + xbm * yrm - (xrm * ybm + xbm * ygm + xgm * yrm)

taur = Kray * Hscl * rho0' should include polarization reduction
tauR0 = taur / Sin(1.486 * drc)'  Sin(1.619 * drc)

'**********************************************
' CHOOSE TROPOSPHERIC AND VOLCANIC AEROSOLS '**
HVul = 2000#                                '**
hV0 = hscroll9.Value * 1000' = 25000#       '**
tauVul = hscroll10.Value / 100              '**
tauVul = tauVul / .0891 / 2                 '**
kLV = -1#                                   '**
beta = hscroll8.Value / 10' Turbidity >=1   '**
kLA = 1#                                    '**
'**********************************************

'****************************************************
' CHOOSE OZONE LAYER HEIGHTS AND THICKNESS        '**
htrop = 15000# 'Tropopause                        '**
hZmx = 26000#                                     '**
nZtrop = .4 ' Tropospheric  Concentration         '**
nZmx = 4    ' Stratospheric Peak                  '**
HZ = 3000   'Ozone Scale Height                   '**
If OZONE = 1 Then sigmZ = .00000051 Else sigmZ = 0'**
Dobson = hscroll11.Value                          '**
sigmZ = sigmZ * Dobson / 300                      '**
kZ = 161.4                                        '**
lmdZ = .59                                        '**
'****************************************************

For il = 1 To 61
lmdR(il) = lmd(il) ^ 4
lmdA(il) = lmd(il) ^ kLA
lmdV(il) = lmd(il) ^ kLV
sgmz(il) = sigmZ * Exp(-kZ * (lmd(il) - lmdZ) ^ 2)
Next il

Select Case kLV
Case -1
Kvul = 25.7
aVul = 3.3: cV = 2 * (1 + aVul ^ 2) / (1 + Exp(-aVul * pi))
Case 0
Kvul = 14#
Case 1
Kvul = 7.57
Case 2
Kvul = 3.96
Case 3
Kvul = 2.1
End Select
bVul = 1 - Exp(-aVul * pi / 2)
Kvul = Kvul * Kray * Hscl / HVul * tauVul

Select Case kLA
Case -1
KAer = 25.7
Case 0
KAer = 14#
Case 1
KAer = 7.57
aAer = 1.5: cA = 2 * (1 + aAer ^ 2) / (1 + Exp(-aAer * pi))
Case 2
KAer = 3.96
Case 3
KAer = 2.1
End Select
bAer = 1 - Exp(-aAer * pi / 2)
KAer = KAer * Kray * Hscl / HAer * (beta - 1)
TauA = HAer * KAer * rho0
TauA0 = TauA / Sin(.821 * drc)' assumes kLa = 1   old value Sin(.9155 * drc)

'******************************************
' START REFRACTION OF VIEWER BEAMS
'******************************************
If REFRACT = 1 Then
 For iphi = 1 To iphimax
 dxp = 1000 * Sgn(Cos(phi(iphi)))
 If phi(iphi) < 1.5 Or phi(iphi) > pi - 1.5 Then
  zf = Re0
  phif = 0
  ix = 0
  hf = h0
   Do Until hf > htop Or hf < 0
    ix = ix + 1
    zf = zf + dxp * Tan(phi(iphi) + phif)
    xf = ix * dxp
    hf = Sqr(zf ^ 2 + xf ^ 2) - Re
    p = p0 * Exp(-hf / Hscl)
    rfi = n0 * p / p0 / Hscl
    phif = phif - dxp * Cos(phi(iphi)) * rfi
   Loop
  xp0 = ix * dxp
  dxf(iphi) = xp0 / ixmax
  zf = Re0
  hfp(ixmax, phi(iphi)) = zf - Re
  phif = 0
  pold = p0
   For ix = 1 To ixmax - 1
    zf = zf + dxf(iphi) * Tan(phi(iphi) + phif)
    hfp(ixmax - ix, iphi) = Sqr(zf ^ 2 + (ix * dxf(iphi)) ^ 2) - Re
    p = p0 * Exp(-hfp(ixmax - ix, iphi) / Hscl)
    rfi = n0 * (pold + p) / Hpscl
    phif = phif - dxf(iphi) * Cos(phi(iphi)) * rfi
    pold = p
   Next ix
 End If
 
 Next iphi
End If
'******************************************
' END REFRACTION OF VIEWER BEAMS
' START SUN ANGLE LOOP  tht0
'******************************************
For ith = 1 To ithmax
tht0 = (tht00 - dth * ith) * drc
x0 = Re0 * Sin(tht0)
z0 = Re0 * Cos(tht0) - Re
dxcld = (hcld - h0) / Tan(phicld) * Cos(tht0) ' FLAT EARTH APPX 22 May 2004
' Determine Cloud Albedo and Transmission
albC = alb(1) + dalb * (1.56 - Abs(tht0)) ^ pwr
trnC = 1 - albC
albE = .2 * albC'Assert Cloud Edge is Optically Thin fraction of albC
trnE = 1 - albE
'******************************************
'START VIEWER BEAM ANGLE LOOP   phi
'******************************************
For iphi = 1 To iphimax
psi = Abs(tht0 - phi(iphi))
If psi > pi Then psi = 2 * pi - psi
If (phi(iphi) < 1.5 Or phi(iphi) > pi - 1.5) And REFRACT = 1 Then
  dtht = Atn(dxf(iphi) * ixmax / Sqr(Retop ^ 2 - (ixmax * dxf(iphi)) ^ 2))
  Else
  dtht = Sgn(pi / 2 - phi(iphi)) * Sqr(Tan(phi(iphi)) ^ 2 + 2 * (htop - h0) / Re) - Tan(phi(iphi))
  If Abs(phi(iphi) - pi / 2) < pi / 4 Then dtht = (htop - h0) / Tan(phi(iphi)) / (Re + h0)
End If
xp0 = Retop * Sin(tht0 + dtht)
zp0 = Retop * Cos(tht0 + dtht) - Re
' x value where ray intersects earth's shadow (zp=0)
dxp = (x0 - xp0) / ixmax
If SHADOW = 1 Then
 If tht0 > phicld Then
   xpen = (hcld - h0) * (1 / Tan(phicld) - (1 - phi(iphi) / phicld) / Tan(tht0))
   Else
   xpen = 0
 End If
End If
If ECLIPSE = 1 Or SHADOW = 1 Then
  ixpen = ixmax - Int(Abs(xpen / dxp))
  If ixpen > ixmax Then ixpen = ixmax
End If
dsp = Abs(dxp / Cos(phi(iphi) - tht0 - dtht))'** REFRACT EQ**
If iphi = iphicld Then ixcld = ixmax - Int(-dxcld / dxp) - 1
'*********************************************
' LOOP FOR SKYLIGHT RAYS
'*********************************************
For ix = 1 To ixmax
xp = xp0 + ix * dxp
If (phi(iphi) < 1.5 Or phi(iphi) > pi - 1.5) And REFRACT = 1 Then
hp = hfp(ix, iphi)
zp = Sqr((Re + hp) ^ 2 - xp ^ 2) - Re
Else
zp = zp0 + ix * dxp * (zp0 - z0) / (xp0 - x0)
End If
thtp = Atn(xp / (zp + Re))' scattering point angle
If REFRACT = 1 And xp < 0 Then
  zss = 328 * Trat ^ 2.9 + 2312 * Trat ^ 1.77 * (xp / 100000#) + 3244 * Exp(2.12 * (Trat - 1)) * tanh(-xp / 248000#)
  Else zss = 0'***1
End If
If zp + Re < 0 Then thtp = pi - thtp
hp = (zp + Re) / Cos(thtp) - Re
If hp < htrop Then nZ = nZtrop Else nZ = nZmx * Exp(-((hp - hZmx) ^ 2 / (hp * HZ)))
' SUNBEAM BRIGHTNESS -
' correct for effective thta near horizon
If zp > zss Or thtp > thtss Then
    If REFRACT = 0 Then
    thtcor = 1.619 * drc * Exp(-((Abs(thtp) / (6# * drc)) ^ .83)) / Trat ^ .7' no refraction
    Else
    thtcor = 1.486 * drc * Exp(-((Abs(thtp) / (5.8 * drc)) ^ .77)) / Trat ^ .7' refraction
    End If
    sgRAY = taur * Exp(-hp / Hscl) / (Sin(Abs(thtp) + thtcor))
    If REFRACT = 0 Then
    thtcor = .9055 * drc * Exp(-((Abs(thtp) / (3.4 * drc)) ^ .8))' no refraction
    Else
    thtcor = .821 * drc * Exp(-((Abs(thtp) / (3.2 * drc)) ^ .77))' refraction
    
    End If
    sgAER = TauA * Exp(-hp / HAer) / (Sin(Abs(thtp) + thtcor))
    If thtp < 0 Then
      sgRAY = 2 * tauR0 * Exp(-zp / Hscl) - sgRAY
      sgAER = 2 * TauA0 * Exp(-zp / HAer) - sgAER
    End If
' Efective Optical Thickness of Sunbeam due to Refraction
    If REFRACT = 1 Then
     If xp < 200000# And zp < 32000 Then
      zReff = (zp - (10500 - 100 * T0 * Exp(xp / (2650 * T0)))) / (66 * T0 - .0004 * xp)
      zAeff = (zp - (10500 - 100 * T0 * Exp(xp / (2650 * T0)))) / (19000 * Exp(xp / 1800000#) - 1000)
      tauRrat = .5 * (1 + tanh(zReff))
      tauArat = .5 * (1 + tanh(zAeff))
      sgRAY = sgRAY * tauRrat
      sgAER = sgAER * tauArat
     End If
    End If
' Calculate tau of sunbeam due to Volcanic Aerosols and Ozone
    If zp > htop Then zp = htop' done 29 April 2005 may cause errors!!!
    xtop = Sqr(htop ^ 2 + 2 * (htop - zp) * Re - zp ^ 2)
    dx2 = (xtop - xp) / (ixmax / 2)
    cc = zp ^ 2 + 2 * zp * Re
   tauVSun = 0
   TauZsun = 0
    For ix2 = 1 To ixmax / 2
     xp2 = xp - dx2 / 2 + ix2 * dx2
     hp2 = -Re + Sqr(Re ^ 2 + (cc + xp2 ^ 2))
     tauVSun = tauVSun + Exp(-Abs(hp2 - hV0) / HVul)
     If hp2 < htrop Then nZ2 = nZtrop Else nZ2 = nZmx * Exp(-((hp2 - hZmx) ^ 2 / (hp2 * HZ)))
     TauZsun = TauZsun + nZ2
    Next ix2
    TauZsun = TauZsun * dx2
    tauVSun = tauVSun * Kvul * dx2
End If
'*** scattering integrated loss vs focused gain (in)
dtauVc = Kvul * Exp(-Abs(hp - hV0) / HVul)
skyRc = Kray * Exp(-hp / Hscl) * dsp
skyAc = KAer * Exp(-hp / HAer) * dsp
dtauVcin = cV * (bVul * Exp(-aVul * psi) + (1 - bVul) * Exp(-aVul * (pi - psi)))
skyRcin = 3 / 4 * (1 + Cos(psi) ^ 2)
skyAcin = cA * (bAer * Exp(-aAer * psi) + (1 - bAer) * Exp(-aAer * (pi - psi)))

'***Wavelength loop to find Optical Thickness TA
For il = 1 To 61
dtauV = dtauVc / lmdV(il)
dtauZ = sgmz(il) * nZ
skyR = skyRc / lmdR(il)
skyA = skyAc / lmdA(il)
skyV = dtauV * dsp
skyZ = dtauZ * dsp
dtauout = skyR + skyA + skyV + skyZ
dtauin = skyR * skyRcin + skyA * skyAcin + skyV * dtauVcin
  If zp > zss Or thtp > thtss Then
   tau = tauVSun / lmdV(il) + TauZsun * sgmz(il) + sgRAY / lmdR(il) + sgAER / lmdA(il)
   F(il) = F0(il) * Exp(-tau)
    If ECLIPSE = 1 Then
      If ix > ixpen Then
      F(il) = .001 * F(il)
      Else
         If (xp - xpen) * 2 / (pi * xpen) > 1 Then
           eclfac = 1
           Else eclfac = (xp - xpen) * 2 / (pi * xpen)
         End If
      F(il) = F(il) * eclfac'((ixpen + 1 - ix) / ixpen)
      End If
    End If

     If ix = ixcld And iphi = iphicld Then
      If phi(iphi) > pi / 2 Or tht0 < 0 Then
        cldil = albC
        edgil = albE
        Else
        cldil = trnC
        edgil = trnE
      End If
      Fcldchrm(il) = F(il) * Exp(.3 * tau * MULT) * cldil * 2' ASSUME FACTOR 2 FOCUSING SINCE EMERGING LIGHT COVERS 1/2 THE SKY
      Fedgchrm(il) = F(il) * Exp(.3 * tau * MULT) * edgil * 2' SINGLE PARTICLE FOCUSING REMOVED 30 OCT 2004
      IrrCld = IrrCld + Fcldchrm(il): IrrEdg = IrrEdg + Fedgchrm(il): IrrSun = IrrSun + F(il)
    End If
    If SHADOW = 1 Then
      If ix > ixpen Then F(il) = .01 * F(il)
    End If

' NEW SKYLIGHT = PENETRATING SKYLIGHT + SCATTERED SUNLIGHT
  Fsct = dtauin * F(il) * Exp(.3 * tau * MULT)
  Else
  Fsct = 0
  dtauout = dtauout - .3 * MULT * (skyR + skyA + skyV)' Mult scatter in shadow
  End If
Fskychrm(il, iphi) = Fsct + Fskychrm(il, iphi) * (1 - dtauout)
    If ix = ixcld And iphi = iphicld Then
      Fcldchrm(il) = Fcldchrm(il) + Fskychrm(il, iphi) * trnpC' need focusing (P) SHOULD THIS BE 2?
      Fedgchrm(il) = Fedgchrm(il) + Fskychrm(il, iphi) * trnpE' need focusing (P) ASKED 30 OCT 2004
    End If
  If ix >= ixcld And iphi = iphicld Then
    Fcldchrm(il) = Fsct + Fcldchrm(il) * (1 - dtauout)
    Fedgchrm(il) = Fsct + Fedgchrm(il) * (1 - dtauout)
  End If
Next il
    If ix = ixcld And iphi = iphicld Then
      zcloud = zp
      thtcld = thtp
    End If
Next ix
'********* END of Sky Beam Loop **********

'Equal Energy Spectrum Rendition
For il = 1 To 61
If Fskychrm(il, iphi) < 0 Then Fskychrm(il, iphi) = 0
xt = xt + Fskychrm(il, iphi) * xbar(il)
yt = yt + Fskychrm(il, iphi) * ybar(il)
zt = zt + Fskychrm(il, iphi) * zbar(il)
irrad(iphi) = irrad(iphi) + xt + yt + zt
Next il
currentx = 0: currenty = 7800
If irrad(iphi) > 0 Then
xchrm(iphi) = xt / (xt + yt + zt)
ychrm(iphi) = yt / (xt + yt + zt)
Else
xchrm(iphi) = 0#
ychrm(iphi) = 0#
End If
If irrad(iphi) > Fmax Then Fmax = irrad(iphi)
xt = 0: yt = 0: zt = 0
Next iphi
'***** SUNBEAM
If z0 >= zss Or tht0 >= thtss Then
 For il = 1 To 61
  xsun = xsun + F(il) * xbar(il)
  ysun = ysun + F(il) * ybar(il)
  zsun = zsun + F(il) * zbar(il)
  Fsun = Fsun + xsun + ysun + zsun
 Next il
 End If
 If z0 >= zss Or tht0 >= thtss Then
  xsunchrm = xsun / (xsun + ysun + zsun)
  ysunchrm = ysun / (xsun + ysun + zsun)
 Else
  xsunchrm = 0
  ysunchrm = 0
 End If
xsun = 0: ysun = 0: zsun = 0
'***** CLOUDBEAMS
      If CLOUD = 1 Then
 For il = 1 To 61
  xcld = xcld + Fcldchrm(il) * xbar(il)
  ycld = ycld + Fcldchrm(il) * ybar(il)
  zcld = zcld + Fcldchrm(il) * zbar(il)
  Fcld = Fcld + xcld + ycld + zcld
  xedg = xedg + Fedgchrm(il) * xbar(il)
  yedg = yedg + Fedgchrm(il) * ybar(il)
  zedg = zedg + Fedgchrm(il) * zbar(il)
  Fedg = Fedg + xedg + yedg + zedg
  Next il
  xcldchrm = xcld / (xcld + ycld + zcld)
  ycldchrm = ycld / (xcld + ycld + zcld)
 xcld = 0: ycld = 0: zcld = 0
  xedgchrm = xedg / (xedg + yedg + zedg)
  yedgchrm = yedg / (xedg + yedg + zedg)
 xedg = 0: yedg = 0: zedg = 0
 If Fedg > Fcld Then FMax2 = Fedg Else FMax2 = Fcld
 If Fmax > FMax2 Then FMax2 = Fmax
 End If
 '***** CHROMATICITY DIAGRM
If CHRM = 1 Then
Line (-1000, 0)-(14000, 8000), RGB(255, 255, 255), BF
  Line (0, 0)-(0, mf)
  Line (0, 0)-(mfx * mf, 0)
  Line (mfx * mf / 3#, mf / 3#)-(mfx * mf * xrm, mf * yrm)
  Line (mfx * mf / 3#, mf / 3#)-(mfx * mf * xgm, mf * ygm)
  Line (mfx * mf / 3#, mf / 3#)-(mfx * mf * xbm, mf * ybm)
  For iphi = 1 To 10
  Line (mfx * mf * .6, mf * (.55 + .03 * iphi))-(mfx * 1.4 * mf, mf * (.55 + .03 * iphi))
  Next iphi
  Line (mfx * mf * .6, mf * .55)-(mfx * 1.4 * mf, mf * .85), , B
  Line (mfx * mf * .6, mf * .7)-(mfx * 1.4 * mf, mf * .7)
  Line (mfx * mf * .8, mf * .55)-(mfx * mf * .8, mf * .85)
  Line (mfx * mf * 1#, mf * .55)-(mfx * mf * 1#, mf * .85)
  Line (mfx * mf * 1.2, mf * .55)-(mfx * mf * 1.2, mf * .85)
  currentx = mfx * mf * .93: currenty = mf * .88
  fontname = "symbol"
  Print "q";
  fontname = "Ms Sans Serif"
  Print "(sun) ="; CInt(10 * tht0 / drc) / 10
  currentx = mfx * mf * .59: currenty = mf * .54
  Print "0                       45                       90                      135                     180"
  currentx = mfx * mf * .72
  Print "Observer Elevation Angle (From Horizon Facing Sun)"
  currentx = mfx * mf * .56: currenty = mf * .56
  Print "0.0"
  currentx = mfx * mf * .56: currenty = mf * .71
  Print "0.5"
  currentx = mfx * mf * .56: currenty = mf * .86
  Print "1.0"
  currentx = mfx * mf * .46: currenty = mf * .8
  Print "Brightness"
  currentx = mfx * mf * .46: currenty = mf * .65
  Print "Color Purity"
  currentx = mfx * mf * .9: currenty = -mf * .01
  Print "X"
  currentx = -300: currenty = mf * .95
  Print "Y"
  
  For il = 1 To 61
    If il > 1 Then Line (mfx * mf * xv(il - 1), mf * yv(il - 1))-(mfx * mf * xv(il), mf * yv(il))
    If il = 61 Then Line (mfx * mf * xv(1), mf * yv(1))-(mfx * mf * xv(61), mf * yv(61))
    thv(il) = Atn((yv(il) - 1# / 3#) / (xv(il) - 1# / 3#))
    If xv(il) - 1# / 3# < 0 Then thv(il) = thv(il) + pi
    If thv(il) < 0 Then thv(il) = thv(il) + 2# * pi
    Next il
 End If
'CALCULATING SKYCOLOR and PURITY
For iphi = 1 To iphimax
Fnorm = irrad(iphi) / Fmax
xrel = xchrm(iphi) - 1# / 3#
yrel = ychrm(iphi) - 1# / 3#
  If CHRM = 1 Then
   pur = Sqr((xchrm(iphi) - 1# / 3#) ^ 2 + (ychrm(iphi) - 1# / 3#) ^ 2)
   thcie = Atn(yrel / xrel)
   If xrel < 0 Then thcie = thcie + pi
   If thcie < 0 Then thcie = thcie + 2 * pi
     If thcie < thv(61) And thcie > thv(1) Then
      xpur = (bp + 1# / 3# * (1# - Tan(thcie))) / (mp - Tan(thcie))
      ypur = mp * xpur + bp
      rp = Sqr((xpur - 1# / 3#) ^ 2 + (ypur - 1# / 3#) ^ 2)
      ElseIf thcie > thv(61) Then
      il = 61
       Do Until thcie < thv(il) Or il = 2
       il = il - 1
       Loop
      rp = Sqr(((xv(il) + xv(il + 1)) / 2# - 1# / 3#) ^ 2 + ((yv(il) + yv(il + 1)) / 2# - 1# / 3#) ^ 2)
      Else
      il = 1
       Do Until thcie > thv(il) Or il = 61
       il = il + 1
       Loop
      rp = Sqr(((xv(il) + xv(il - 1)) / 2# - 1# / 3#) ^ 2 + ((yv(il) + yv(il - 1)) / 2# - 1# / 3#) ^ 2)
     End If
    purpct = pur / rp
   End If
xcolor = Sgn(xrel) * Abs(xrel) ^ (1 + .3 * (1 - CHRM)) + xneut'xrel / 2 + xneut
ycolor = Sgn(yrel) * Abs(yrel) ^ (1 + .3 * (1 - CHRM)) + yneut  'yrel / 2 + yneut
detR = xcolor * ygm + xgm * ybm + xbm * ycolor - (xcolor * ybm + xbm * ygm + xgm * ycolor)
detG = xrm * ycolor + xcolor * ybm + xbm * yrm - (xrm * ybm + xbm * ycolor + xcolor * yrm)
detB = xrm * ygm + xgm * ycolor + xcolor * yrm - (xrm * ycolor + xcolor * ygm + xgm * yrm)
r = detR / detx: g = detG / detx: B = detB / detx
If r > g Then huemax = r Else huemax = g
If B > huemax Then huemax = B
If r < 0 Then r = 0
If g < 0 Then g = 0
If B < 0 Then B = 0
If r = 0 And g = 0 Then B = 0
' logarithmic brightness scale (preserving color)
r = Int(255 * r / huemax)
g = Int(255 * g / huemax)
B = Int(255 * B / huemax)
r = CInt(r * Fnorm ^ eye)
g = CInt(g * Fnorm ^ eye)
B = CInt(B * Fnorm ^ eye)
 If CHRM = 1 Then
  fillstyle = 0
  fillcolor = RGB(r, g, B)
  If Fnorm > 0 Then
  Circle (mfx * mf * xchrm(iphi), mf * ychrm(iphi)), 20 + 2 * (iphimax - iphi)
  Circle (mfx * mf * (.6 + .8 * phi(iphi) / pi), mf * (.55 + .3 * purpct)), 40 + 2 * (iphimax - iphi)
  End If
  fillcolor = RGB(0, 0, 0)
  Circle (mfx * mf * (.6 + .8 * phi(iphi) / pi), mf * (.55 + .3 * Fnorm ^ .5)), 30
  fillstyle = 1
 End If
 If CHRM = 0 Then
'  sun color here Circle (2500, 3000), 1000, RGB(R, G, B)
   Line (xfrm * (ith - 1), czdraw * phi(iphi))-(xfrm * ith, czdraw * phi(iphi - 1)), RGB(r, g, B), BF
   If iphi = 1 Then Line (xfrm * (ith - 1), czdraw * phi(1))-(xfrm * ith, czdraw * (phi(1) - .01)), RGB(r, g, B), BF
   End If
huemax = 0#
For il = 1 To 61
Fskychrm(il, iphi) = 0
Next il
 If CHRM = 0 Then
   currentx = 0
   'currenty = 8000 - 175 * iphi
   'Print CInt(phi(iphi) / drc * 100) / 100, CLng(irrad(iphi)), r, g, B
      IrrCld = 0: IrrEdg = 0: IrrSun = 0
End If
irrad(iphi) = 0
Next iphi
If CHRM = 1 Then
currentx = 8800
currenty = 6900
End If
'fontname = "Symbol"
'Print "q";
'fontname = "Ms Sans Serif"
'Print "(sun) = "; CInt(tht0 / drc)
'***** SUN, CLOUD, AND CLOUD EDGE COLOR
For i = 1 To 3
If i = 1 Then
   If z0 >= zss Or tht0 >= thtss Then
    xrel = xsunchrm - 1# / 3#
    yrel = ysunchrm - 1# / 3#
   End If
  ElseIf i = 2 Then
   xrel = xedgchrm - 1# / 3#
   yrel = yedgchrm - 1# / 3#
  Else
   xrel = xcldchrm - 1# / 3#
   yrel = ycldchrm - 1# / 3#
End If
xcolor = Sgn(xrel) * Abs(xrel) ^ (1 + .3 * (1 - CHRM)) + xneut
ycolor = Sgn(yrel) * Abs(yrel) ^ (1 + .3 * (1 - CHRM)) + yneut
detR = xcolor * ygm + xgm * ybm + xbm * ycolor - (xcolor * ybm + xbm * ygm + xgm * ycolor)
detG = xrm * ycolor + xcolor * ybm + xbm * yrm - (xrm * ybm + xbm * ycolor + xcolor * yrm)
detB = xrm * ygm + xgm * ycolor + xcolor * yrm - (xrm * ycolor + xcolor * ygm + xgm * yrm)
r = detR / detx: g = detG / detx: B = detB / detx
If r > g Then huemax = r Else huemax = g
If B > huemax Then huemax = B
If r < 0 Then r = 0
If g < 0 Then g = 0
If B < 0 Then B = 0
r = Int(255 * r / huemax)
g = Int(255 * g / huemax)
B = Int(255 * B / huemax)
If i = 1 And (z0 >= zss Or tht0 >= thtss) Then
   If CHRM = 1 Then
     fillstyle = 0
     fillcolor = RGB(r, g, B)
     'Circle (mfx * mf * xsunchrm, mf * ysunchrm), 150
     Line (mfx * mf * xsunchrm - 100, mf * ysunchrm - 100)-(mfx * mf * xsunchrm + 100, mf * ysunchrm + 100), RGB(r, g, B), BF
     fillstyle = 1
   End If
   If CHRM = 0 Then
   fillstyle = 0
   fillcolor = RGB(r, g, B)
   Circle (xfrm * (ith - .5), czdraw * tht0), rsun, RGB(r, g, B)
   End If
   ElseIf i = 2 And CLOUD = 1 Then
     r = CInt(r * (Fedg / FMax2) ^ eye)
     g = CInt(g * (Fedg / FMax2) ^ eye)
     B = CInt(B * (Fedg / FMax2) ^ eye)
currentx = 0: currenty = 8000 - 150 * (iphimax + 6)
   If CLOUD = 1 And CHRM = 1 Then
     fillstyle = 0
     fillcolor = RGB(r, g, B)
     Line (mfx * mf * xedgchrm - 150, mf * yedgchrm - 150)-(mfx * mf * xedgchrm + 150, mf * yedgchrm + 150), , B
     fillstyle = 1
   End If
 If CLOUD = 1 And CHRM = 0 Then
   drawwidth = 1
   fillcolor = RGB(r, g, B)
   fillstyle = 0
  For j = 0 To 2
   For ix = 4 To 4 - j Step -1
     Circle (xfrm * (ith - 1.4) + 1.5 * rcld * ix + rcld / 2 * j, .75 * rcld * j + czdraw * (.002 + phicld)), rcld, RGB(r, g, B), , , .33
   Next ix
  Next j
 End If
   ElseIf CLOUD = 1 Then
      r = CInt(r * (Fcld / FMax2) ^ eye)
      g = CInt(g * (Fcld / FMax2) ^ eye)
      B = CInt(B * (Fcld / FMax2) ^ eye)
      currentx = 0: currenty = 8000 - 150 * (iphimax + 7)
   If CLOUD = 1 And CHRM = 1 Then
     fillstyle = 0
     fillcolor = RGB(r, g, B)
     Line (mfx * mf * xcldchrm - 100, mf * ycldchrm - 100)-(mfx * mf * xcldchrm + 100, mf * ycldchrm + 100), , B
     fillstyle = 1
   End If
   If CLOUD = 1 And CHRM = 0 Then
    fillcolor = RGB(r, g, B)
   For j = 0 To 2
     For ix = 4 To 4 - j Step -1
     Circle (xfrm * (ith - 1.4) + 1.5 * rcld * ix + rcld / 2 * j, -40 + .75 * rcld * j + czdraw * phicld), rcld, RGB(r, g, B), , , .1
     Circle (xfrm * (ith - 1.4) + 1.5 * rcld * ix + rcld / 2 * j, -70 + .75 * rcld * j + czdraw * phicld), rcld, RGB(r, g, B), , , .15
                                                                                 
   Next ix
    Next j
  End If
End If
Next i
'*****
drawwidth = 1
fillstyle = 1
If CHRM = 0 Then
  Line (xfrm * (ith - 1), czdraw * (phi(1) - .01))-(xfrm * ith, czdraw * phi(iphimax)), RGB(0, 0, 0), B
  currentx = xfrm * (ith - .9): currenty = czdraw * phi(iphimax) + 300
  fontname = "symbol"
  Print "q";
  fontname = "Ms Sans Serif"
  Print "(sun) = "; CInt(10 * tht0 / drc) / 10
End If
' End of Equal Energy Color Analysis
Fcld = 0: Fedg = 0
currentx = 0
Next ith
Close #1
End Sub

Sub HScroll1_Change ()
Dim tht00 As Single
Cls
currentx = 3000
currenty = 1000
tht00 = hscroll1.Value
Print tht00;
Print " Initial Sun Elev Angle"
End Sub

Sub HScroll10_Change ()
Dim tauVul As Single
Cls
currentx = 3000
currenty = 1000
tauVul = hscroll10.Value / 100
Print tauVul;
Print " Volcanic Aerosol Optical Thickness"
End Sub

Sub HScroll11_Change ()
Dim Dobson As Single
Cls
currentx = 3000
currenty = 1000
Dobson = hscroll11.Value
Print Dobson;
Print " Ozone Depth (Dobson Units)"
End Sub

Sub HScroll12_Change ()
Dim SKYFRAC As Single
Cls
currentx = 3000
currenty = 1000
SKYFRAC = hscroll12.Value / 20
Print SKYFRAC;
Print " Fraction of Sky Seen (0.5 = Horizon to Zenith)"
End Sub

Sub HScroll13_Change ()
Dim T0 As Single
Cls
currentx = 3000
currenty = 1000
T0 = hscroll13.Value
Print T0;
Print " Atmospheric Temperature"
End Sub

Sub HScroll2_Change ()
Dim dtht As Single
Cls
currentx = 3000
currenty = 1000
dtht = hscroll2.Value / 10
Print dtht;
Print " Increment Sun Elevation Angle"
End Sub

Sub HScroll3_Change ()
Dim h0 As Single
Cls
currentx = 3000
currenty = 1000
h0 = hscroll3.Value * 100
Print h0;
Print " Observer height (m)"
End Sub

Sub HScroll4_Change ()
Dim hcld As Single
Cls
currentx = 3000
currenty = 1000
hcld = hscroll4.Value * 100
Print hcld;
Print " Cloud height (m) above observer"
End Sub

Sub HScroll5_Change ()
Dim phicld As Single
Cls
currentx = 3000
currenty = 1000
phicld = hscroll5.Value
Print phicld;
Print " Cloud Elev Angle"
End Sub

Sub HScroll6_Change ()
Dim taucld As Single
Cls
currentx = 3000
currenty = 1000
taucld = hscroll6.Value
Print taucld;
Print " Cloud optical thickness "
End Sub

Sub HScroll7_Change ()
Dim HAer As Single
Cls
currentx = 3000
currenty = 1000
HAer = hscroll7.Value * 100
Print HAer;
Print " Aerosol Scale Height (m)"
End Sub

Sub HScroll8_Change ()
Dim beta As Single
Cls
currentx = 3000
currenty = 1000
beta = hscroll8.Value / 10
Print beta;
Print " Atmospheric turbidity (Pure Air = 1)"
End Sub

Sub HScroll9_Change ()
Dim hV0 As Single
Cls
currentx = 3000
currenty = 1000
hV0 = hscroll9.Value * 1000
Print hV0;
Print " Center Height Volcanic Aerosols (m)"
End Sub

Sub Label1_Click ()
End Sub

Sub Option1_Click ()
End Sub

Sub Check2_Click ()
End Sub