Solar Radiation Model for QuickBasic

R. H. B. Exell, A program in BASIC for calculating solar radiation in tropical climates on small computers. Renewable Energy Review Journal, Vol. 8, No. 2, December 1986, Regional Energy Resources Information Center, Asian Institute of Technology, Bangkok.

Trigonometric Calculations

FUNCTION Arccos
FUNCTION Arcsin
FUNCTION Arctan

Finding the Sun's Position

FUNCTION Date
SUB Declination
SUB HourAngle
SUB SolarPosition

Finding Solar Radiation

SUB SolarIrradiance
SUB TiltedSurfaceIrradiance

Daily Solar Radiation in Changing Weather

SUB MeanDailyIrradiation
SUB DailyClearnessDistribution
SUB RandomDailyClearness

Hourly Solar Radiation

SUB HourlySolarPosition
SUB RandomHourlyClearness
SUB HourlySolarIrradiance
SUB HourlyTiltedSurfaceIrradiance

REM  Solar Radiation Model for QBasic.
'    Copyright 1997, R. H. B. Exell
'    King Mongkut's Institute of Technology Thonburi

REM  This model is for calculating solar radiation
'    falling on objects at the surface of the Earth
'    in the tropics.  Write your own main program
'    to do the calculations you want.  Use the module-level
'    variables to transfer data to and from SUB procedures.

DIM SHARED Phi AS SINGLE, Lambda AS SINGLE, Zone AS SINGLE
REM  Phi = Latitude (deg N).
'    Lambda = Longitude (deg E).
'    Zone = Time zone (hours ahead of GMT).

DIM SHARED Nd AS SINGLE, ZMT AS SINGLE, Eqt AS SINGLE
REM  Nd = Day in the year (1 Jan = 1, ..., 31 Dec = 365).
'    ZMT = Zone Mean Time (hours from midnight).
'    Eqt = Equation of time (min).

DIM SHARED Delta AS SINGLE, Omega AS SINGLE, Psi AS SINGLE, Zeta AS SINGLE
REM  Delta = Declination of the sun (deg N).
'    Omega = Hour angle of the sun (deg W).
'    Psi = Azimuth of the sun (deg W of S).
'    Zeta = Zenith angle of the sun (deg).

DIM SHARED Ksky AS SINGLE, Ig AS SINGLE, Ib AS SINGLE, Id AS SINGLE
REM  Ksky = Clearness of the sky (clear sky = 1).
'    Ig = Global solar irradiance (kW/m2).
'    Ib = Beam solar irradiance (kW/m2).
'    Id = Diffuse solar irradiance (kW/m2).

DIM SHARED Alpha AS SINGLE, Beta AS SINGLE, CosTheta AS SINGLE, It AS SINGLE
REM  Alpha = Azimuth of a tilted surface (deg W of S).
'    Beta = Tilt angle of a surface (deg).
'    CosTheta = Cosine of the incidence angle.
'    It = Solar irradiance on a tilted surface (kW/m2).

DIM SHARED Hmax AS SINGLE, Hmean AS SINGLE, Kmean AS SINGLE, Kd AS SINGLE
REM  Hmax = Maximum daily global solar irradiation (MJ/m2).
'    Hmean = Mean daily global solar irradiation (MJ/m2).
'    Kmean = Mean daily clearness of the sky (clear sky = 1).
'    Kd = Daily clearness of the sky (clear sky = 1).

DIM SHARED Hmonth(13) AS SINGLE, PKd(10) AS SINGLE
REM  Hmonth() = Monthly mean daily irradiation (MJ/m2);
'         Hmonth(1) for Jan, ..., Hmonth(12) for Dec.
'    PKd() = Cumulative probability distribution of daily
'         clearness values; PKd(k) is the probability
'         that Kd is less than 0.1k.

DIM SHARED Psih(12) AS SINGLE, Zetah(12) AS SINGLE
REM  Psih() = Hourly azimuth of the sun (deg W of S)
'         from 6 a.m. to 6 p.m. apparent solar time.
'    Zetah() = Hourly zenith angle of the sun (deg)
'         from 6 a.m. to 6 p.m. apparent solar time.

DIM SHARED Kh(12) AS SINGLE, Igh(12) AS SINGLE
REM  Kh() = Hourly clearness of the sky
'         from 6 a.m. to 6 p.m. apparent solar time.
'    Igh() = Hourly global solar irradiance (kW/m2)
'         from 6 a.m. to 6 p.m. apparent solar time.

DIM SHARED Ibh(12) AS SINGLE, Idh(12) AS SINGLE
REM  Ibh() = Hourly beam solar irradiance (kW/m2)
'         from 6 a.m. to 6 p.m. apparent solar time.
'    Idh() = Hourly diffuse solar irradiance (kW/m2)
'         from 6 a.m. to 6 p.m. apparent solar time.

DIM SHARED CosThetah(12) AS SINGLE, Ith(12) AS SINGLE
REM  CosThetah() = Hourly cosine of the incidence angle
'         from 6 a.m. to 6 p.m. apparent solar time.
'    Ith() = Hourly solar irradiance on a tilted surface
'         from 6 a.m. to 6 p.m. apparent solar time.

CONST DR = 1.745329E-02
REM  Converts degrees to radians (x! degrees = x! * DR radians).
	  
DECLARE FUNCTION Arcsin! (Sine!)
REM  The angle (deg) with sine equal to "Sine!".

DECLARE FUNCTION Arccos! (Cosine!)
REM  The angle (deg) with cosine equal to "Cosine!".

DECLARE FUNCTION Arctan! (Sine!, Cosine!)
REM  The angle (deg) with sine equal to "Sine!"
'    and cosine equal to "Cosine!".

DECLARE FUNCTION Date% (Day%, Month%)
REM  The day in the year (1 Jan = 1, ..., 31 Dec = 365)
'    on the given date.

DECLARE SUB Declination ()
REM  Finds Delta from Nd.

DECLARE SUB HourAngle ()
REM  Finds Eqt and Omega from Lambda, Zone, Nd, and ZMT.

DECLARE SUB SolarPosition ()
REM  Finds Zeta and Psi from Phi, Delta, and Omega.

DECLARE SUB SolarIrradiance ()
REM  Finds Ig, Ib, and Id from Nd, Zeta, and Ksky.

DECLARE SUB TiltedSurfaceIrradiance ()
REM  Finds CosTheta and It
'    from Alpha, Beta, Zeta, Psi, Ib, and Id.

DECLARE SUB MeanDailyIrradiation ()
REM  Finds Hmean, Hmax, and Kmean
'    from Hmonth(), Phi, and Nd.

DECLARE SUB DailyClearnessDistribution ()
REM  Finds PKd() from Kmean.

DECLARE SUB RandomDailyClearness ()
REM  Finds Kd from PKd().

DECLARE SUB HourlySolarPosition ()
REM  Finds Zetah() and Psih() from Phi and Delta.

DECLARE SUB RandomHourlyClearness ()
REM  Finds Kh() from Kd.

DECLARE SUB HourlySolarIrradiance ()
REM  Finds Igh(), Ibh(), and Idh()
'    from Nd, Kh(), and Zetah().

DECLARE SUB HourlyTiltedSurfaceIrradiance ()
REM  Finds CosThetah() and Ith()
'    from Alpha, Beta, Zetah(), Psih, Ibh(), and Idh().

REM  Insert your program here. Back to where you were

END

FUNCTION Arccos! (Cosine!)
IF ABS(Cosine!) < .999999 THEN
	Arccos! = 90 - 57.29578 * ATN(Cosine! / SQR(1 - Cosine! * Cosine!))
ELSE
	Arccos! = 90 - 90 * SGN(Cosine!)
END IF
END FUNCTION

FUNCTION Arcsin! (Sine!)
IF ABS(Sine!) < .999999 THEN
	Arcsin! = 57.29578 * ATN(Sine! / SQR(1 - Sine! * Sine!))
ELSE
	Arcsin! = 90 * SGN(Sine!)
END IF
END FUNCTION

FUNCTION Arctan! (Sine!, Cosine!)
IF ABS(Cosine!) < .000001 THEN
	Arctan! = 90 * SGN(Sine!)
ELSEIF Cosine! > 0 THEN
	Arctan! = 57.29578 * ATN(Sine! / Cosine!)
ELSEIF Sine! = 0 THEN
	Arctan! = 180
ELSE
	LET x! = 57.29578 * ATN(Sine! / Cosine!)
	Arctan! = x! - 180 * SGN(x!)
END IF
END FUNCTION

SUB DailyClearnessDistribution
LET x! = Kmean
IF x! < .05 THEN x! = .05
IF x! > .95 THEN x! = .95
LET p! = (x! - .05) / .9: q! = 1 - p!: r! = p! * p! * p!: s! = q! * q! * q!
PKd(1) = s! * s! * s!: PKd(2) = 9 * p! * s! * s! * q! * q!
PKd(3) = 36 * p! * p! * s! * s! * q!: PKd(4) = 84 * r! * s! * s!
PKd(5) = 126 * r! * p! * s! * q! * q!: PKd(6) = 126 * r! * p! * p! * s! * q!
PKd(7) = 84 * r! * r! * s!: PKd(8) = 36 * r! * r! * p! * q! * q!
PKd(9) = 9 * r! * r! * p! * p! * q!
PKd(0) = 0
FOR k% = 1 TO 8
	PKd(k% + 1) = PKd(k%) + PKd(k% + 1)
NEXT k%
PKd(10) = 1
END SUB

FUNCTION Date% (Day%, Month%)
IF Month% = 1 THEN Date% = Day%
IF Month% = 2 THEN Date% = 31 + Day%
IF Month% = 3 THEN Date% = 59 + Day%
IF Month% = 4 THEN Date% = 90 + Day%
IF Month% = 5 THEN Date% = 120 + Day%
IF Month% = 6 THEN Date% = 151 + Day%
IF Month% = 7 THEN Date% = 181 + Day%
IF Month% = 8 THEN Date% = 212 + Day%
IF Month% = 9 THEN Date% = 243 + Day%
IF Month% = 10 THEN Date% = 273 + Day%
IF Month% = 11 THEN Date% = 304 + Day%
IF Month% = 12 THEN Date% = 334 + Day%
END FUNCTION

SUB Declination
LET x! = (Nd - 80) * .9863
Delta = .386 + 23.273 * COS((x! - 92) * DR) + .4 * COS((x! + x! - 19.2) * DR)
END SUB

SUB HourAngle
LET x! = (Nd - 80) * .9863
Eqt = .013 + 7.342 * COS((x! - 194.9) * DR) + 9.939 * COS((x! + x! - 93.1) * DR)
Omega = (ZMT - 12 - Zone) * 15 + Lambda + Eqt * .25
END SUB

SUB HourlySolarIrradiance
FOR h% = 0 TO 12
	Zeta = Zetah(h%): Ksky = Kh(h%)
	SolarIrradiance
	Ibh(h%) = Ib: Idh(h%) = Id: Igh(h%) = Ig
NEXT h%
END SUB

SUB HourlySolarPosition
FOR h% = 0 TO 12
	Omega = h% * 15 - 90
	SolarPosition
	Psih(h%) = Psi: Zetah(h%) = Zeta
NEXT h%
END SUB

SUB HourlyTiltedSurfaceIrradiance
FOR h% = 0 TO 12
	Zeta = Zetah(h%): Psi = Psih(h%)
	Ib = Ibh(h%): Id = Idh(h%)
	TiltedSurfaceIrradiance
	CosThetah(h%) = CosTheta: Ith(h%) = It
NEXT h%
END SUB

SUB MeanDailyIrradiation
REM  Calculation of Hmean using monthly data.
Hmonth(0) = Hmonth(12): Hmonth(13) = Hmonth(1)
LET x! = Nd * .0328767 + .5: N% = INT(x!): p! = x! - N%
Hmean = (1 - p!) * Hmonth(N%) + p! * Hmonth(N% + 1)
REM  Calculation of Hmax.
IF ABS(Phi) > 25 THEN
	BEEP
	PRINT "Phi out of range for Hmax and Kmean"
	PRINT "Press any key to continue"
	DO WHILE INKEY$ = ""
	LOOP
END IF
LET x! = (Nd - 80) * .9863
f! = 1 - .0335 * SIN(.9863 * (Nd - 94) * DR)
a! = (-.0041215 * Phi + .00325) * Phi + 27.4486
b! = .321 * Phi
c! = (-.0006025 * Phi + .0024) * Phi + 1.215
Hmax = (a! + b! * COS((x! - 92) * DR) + c! * COS((x! + x! - .12 * Phi - 3.4) * DR)) * f!
REM  Calculation of Kmean.
Kmean = Hmean / Hmax
END SUB

SUB RandomDailyClearness
LET r! = RND: k% = 0
DO WHILE PKd(k% + 1) < r!
	k% = k% + 1
LOOP
Kd = .1 * (k% + (r! - PKd(k%)) / (PKd(k% + 1) - PKd(k%)))
END SUB

SUB RandomHourlyClearness
REM  Calculation of hourly clearness distribution
DIM PKh(10)
LET x! = Kd
IF x! < .05 THEN x! = .05
IF x! > .95 THEN x! = .95
LET p! = (x! - .05) / .9: q! = 1 - p!: r! = p! * p! * p!: s! = q! * q! * q!
PKh(1) = s! * s! * s!: PKh(2) = 9 * p! * s! * s! * q! * q!
PKh(3) = 36 * p! * p! * s! * s! * q!: PKh(4) = 84 * r! * s! * s!
PKh(5) = 126 * r! * p! * s! * q! * q!: PKh(6) = 126 * r! * p! * p! * s! * q!
PKh(7) = 84 * r! * r! * s!: PKh(8) = 36 * r! * r! * p! * q! * q!
PKh(9) = 9 * r! * r! * p! * p! * q!
PKh(0) = 0
FOR k% = 1 TO 8
	PKh(k% + 1) = PKh(k%) + PKh(k% + 1)
NEXT k%
PKh(10) = 1
REM  Generation of random hourly clearness values
DO
	FOR h% = 0 TO 12
		LET r! = RND: k% = 0
		DO WHILE PKh(k% + 1) < r!
			k% = k% + 1
		LOOP
		Kh(h%) = .1 * (k% + (r! - PKh(k%)) / (PKh(k% + 1) - PKh(k%)))
	NEXT h%
	LET s! = .034 * (Kh(1) + Kh(11)) + .066 * (Kh(2) + Kh(10)) + .093 * (Kh(3) + Kh(9)) + .114 * (Kh(4) + Kh(8)) + .127 * (Kh(5) + Kh(7)) + .132 * Kh(6)
LOOP UNTIL ABS(s! - Kd) < .05
END SUB

SUB SolarIrradiance
IF Zeta >= 90 THEN
	Ig = 0: Ib = 0: Id = 0
ELSE
	LET z! = Zeta * Zeta / 8100
	f! = 1 - .0335 * SIN(.9863 * (Nd - 94) * DR)
	Ig = ((((((-4.4631 * z! + 13.0798) * z! - 13.899) * z! + 6.6849) * z! - 1.072) * z! - 1.4354) * z! + 1.1049) * f! * Ksky
	Ib = ((((((-.9217 * z! + 3.7206) * z! - 5.7674) * z! + 3.3705) * z! - 1.1883) * z! - .2001) * z! + .9864) * f! * Ksky
	Id = Ig - Ksky * Ksky * (.733 + Ksky * Ksky * .267) * Ib * COS(Zeta * DR)
	IF Zeta < 87.5 THEN
		Ib = (Ig - Id) / COS(Zeta * DR)
	ELSE
		Ib = 0
	END IF
END IF
END SUB

SUB SolarPosition
LET cp! = COS(Phi * DR): cd! = COS(Delta * DR): co! = COS(Omega * DR)
sp! = SIN(Phi * DR): sd! = SIN(Delta * DR): so! = SIN(Omega * DR)
Zeta = Arccos!(cp! * cd! * co! + sp! * sd!)
IF Zeta = 0 OR Zeta = 180 THEN
	Psi = 0
ELSE
	LET sz! = SIN(Zeta * DR)
	Sine! = cd! * so! / sz!
	Cosine! = (sp! * cd! * co! - cp! * sd!) / sz!
	Psi = Arctan!(Sine!, Cosine!)
END IF
END SUB

SUB TiltedSurfaceIrradiance
CosTheta = COS(Zeta * DR) * COS(Beta * DR) + SIN(Zeta * DR) * SIN(Beta * DR) * COS((Psi - Alpha) * DR)
IF CosTheta <= 0 THEN
	It = .5 * Id * (1 + COS(Beta * DR))
ELSE
	It = Ib * CosTheta + .5 * Id * (1 + COS(Beta * DR))
END IF
END SUB

Home Page   Introduction to QuickBasic

By R. H. B. Exell, 1998. King Mongkut's University of Technology Thonburi.