TheBackShed.com

Forum Index : Microcontroller and PC projects : Sun position and PV output

Author

Message

TassyJim

Guru

Joined: 07/08/2011
Location: Australia
Posts: 6297

Posted: 07:32am 03 Aug 2025

Copy link to clipboard

Print this post

A program to estimate the output of PV arrays at any position. It can be used to compare multiple layouts. Handy when you are deciding where to place the arrays.

You need to enter your lat, long, timezone and the array declination (tilt) and azimuth(bearing).
eff is the estimated efficiency accounting for system losses and cloud cover etc. Using 0.75 agrees with my location.
directions are true, not magnetic.
The program assumes 2 arrays but you could add as many as required.

The program first prints the current days output at 15 minute intervals.
Then it calculates the total output of one day for each week of the year.
The code to chart the year output is commented out.

The noaaSun SUB is a translation of a NOAA spreadsheet.
https://gml.noaa.gov/grad/solcalc/azel.html
It's main purpose is to give the sun elevation and azimuth for any place and time.

The function arrayOutput() returns the maximum output of a PV array.
Before running arrayOutput() you need to call the noaaSun sub which populates the SUN.xxx variables for that time.

I used a simple cos() function to determine the output. In reality, there are many variables that could be taken into account but simple is enough to compare systems.

The ZIP contains the code with line continuation and a second copy with long lines.
Note that one line would be longer than 255 characters so I used intermediate temp variables to allow it to run.

If you use the version with line continuation, you must set OPTION CONTINUATION lines on
If your device or MMB4W doesn't have OPTION CONTINUATION lines, use the sun_angle_LL.bas

'
' OPTION CONTINUATION lines on
OPTION EXPLICIT
OPTION DEFAULT FLOAT
DIM SUN.date,SUN.time,SUN.jd,SUN.jc,SUN.gmLongSun,SUN.gmAnomSun,SUN.eccEarth
DIM SUN.sunEqu,SUN.sunTrueLong,SUN.sunTrueAnom,SUN.suRadV,SUN.sunAppLong
DIM SUN.meanObliptic,SUN.obliqCorr,SUN.sunRtAscen,SUN.sunDec,SUN.y,SUN.equTime
DIM SUN.HAsunrise,SUN.solarNoon,SUN.sunriseTime,SUN.sunsetTime,SUN.sunDuration
DIM SUN.trueSolarTime,SUN.hourAngle,SUN.solarZenithAng,SUN.solarElevAng
DIM SUN.atmosREfr,SUN.solarElevCorr,SUN.SolarAzAngle

DIM lat = -41.12545, long = 145.90739, tz = 10

DIM array1Az = 90, array1Dec = 15 , array1Max = 2.5, array1Total
DIM array2Az = 270, array2Dec = 15, array2Max = 2.5, array2Total
DIM eff = 0.75
DIM p1,p2,t,n,d, tdate
DIM arrayDay(52)
DIM weeks(52)
FOR n = 0 TO 52
weeks(n) = n*6
NEXT n
DIM today$ = DATE$
DIM thistime$ = "00:00:00"

day:
FOR t = 0 TO 1425 STEP 15
thisTime$ = STR$(t\60,2,0,"0")+":"+STR$(t MOD 60,2,0,"0")+":00"
tdate = EPOCH(today$+" "+thistime$) + d*7*86400

noaaSun(tdate, lat, long, tz)
p1 = arrayOutput(array1Az, array1Dec,array1Max)*eff
p2 = arrayOutput(array2Az, array2Dec,array2Max)*eff
array1Total = array1Total + p1/4
array2Total = array2Total + p2/4
'print SUN.date
'print SUN.time
'print SUN.jd
'print SUN.jc
'print SUN.solarElevAng
IF SUN.solarElevAng >= 0 THEN
PRINT thistime$,STR$(p1,2,3),STR$(p2,2,3),STR$(p1+p2,2,3), _
STR$(SUN.solarElevAng,3,1)
ENDIF
NEXT t
PRINT LEFT$(DATETIME$(tdate),10),STR$(array1Total+array2Total,2,3)
PRINT
'end
year:
today$ = "01/01/2025"
thistime$ = "00:00:00"
FOR d = 0 TO 52
array1Total = 0
array2Total = 0
FOR t = 0 TO 1425 STEP 15
thisTime$ = STR$(t\60,2,0,"0")+":"+STR$(t MOD 60,2,0,"0")+":00"
tdate = EPOCH(today$+" "+thistime$) + d*7*86400

noaaSun(tdate, lat, long, tz)
p1 = arrayOutput(array1Az, array1Dec,array1Max)*eff
p2 = arrayOutput(array2Az, array2Dec,array2Max)*eff
array1Total = array1Total + p1/4
array2Total = array2Total + p2/4
NEXT t
PRINT LEFT$(DATETIME$(tdate),10),STR$(array1Total+array2Total,2,3)
arrayDay(d)=array1Total+array2Total
NEXT d
'cls
'math window arrayDay(),200,20,arrayDay()
'line graph weeks(),arrayDay()

SUB noaaSun(tdate, lat, long, tz )
LOCAL tmp1, tmp2
SUN.date = INT(tdate / 86400) + 25569
SUN.time = (tdate MOD 86400) / 86400
SUN.jd = SUN.date+2415018.5+SUN.time-tz/240
SUN.jc = (SUN.jd-2451545)/36525
SUN.gmLongSun = (280.46646+SUN.jc*(36000.76983 + SUN.jc*0.0003032) MOD 360)
SUN.gmAnomSun = 357.52911+SUN.jc*(35999.05029 - 0.0001537*SUN.jc)
SUN.eccEarth = 0.016708634-SUN.jc*(0.000042037+0.0000001267*SUN.jc)
SUN.sunEqu = SIN(RAD(SUN.gmAnomSun))*(1.914602-SUN.jc* _
(0.004817+0.000014*SUN.jc))+SIN(RAD(2*SUN.gmAnomSun))* _
(0.019993-0.000101*SUN.jc)+SIN(RAD(3*SUN.gmAnomSun))*0.000289
SUN.sunTrueLong = SUN.gmLongSun + SUN.sunEqu
SUN.sunTrueAnom = SUN.gmAnomSun + SUN.gmLongSun
SUN.suRadV = (1.000001018*(1-SUN.eccEarth*SUN.eccEarth))/ _
(1+SUN.eccEarth*COS(RAD(SUN.sunTrueAnom)))
SUN.sunAppLong = SUN.sunTrueLong-0.00569-0.00478* _
SIN(RAD(125.04-1934.136*SUN.jc))
SUN.meanObliptic = 23+(26+((21.448-SUN.jc*(46.815+SUN.jc* _
(0.00059-SUN.jc*0.001813))))/60)/60
SUN.obliqCorr = SUN.meanObliptic+0.00256*COS(RAD(125.04-1934.136*SUN.jc))
SUN.sunRtAscen = DEG(ATAN2(COS(RAD(SUN.sunAppLong)),COS(RAD(SUN.obliqCorr)) _
*SIN(RAD(SUN.sunAppLong))))
SUN.sunDec = DEG(ASIN(SIN(RAD(SUN.obliqCorr))*SIN(RAD(SUN.sunAppLong))))
SUN.y = TAN(RAD(SUN.obliqCorr/2))*TAN(RAD(SUN.obliqCorr/2))
tmp1 = SUN.y*SIN(2*RAD(SUN.gmLongSun))-2* SUN.eccEarth* _
SIN(RAD(SUN.gmAnomSun))+4* SUN.eccEarth*SUN.y*SIN(RAD(SUN.gmAnomSun))* _
COS(2*RAD(SUN.gmLongSun))
tmp2 = 0.5*SUN.y*SUN.y*SIN(4*RAD(SUN.gmLongSun))-1.25*SUN.eccEarth* _
SUN.eccEarth*SIN(2*RAD(SUN.gmAnomSun))
SUN.equTime = 4*DEG(tmp1-tmp2)
SUN.HAsunrise = DEG(ACOS(COS(RAD(90.833))/(COS(RAD(lat))* _
COS(RAD(SUN.sunDec)))-TAN(RAD(lat))*TAN(RAD(SUN.sunDec))))
SUN.solarNoon = (720-4*long-SUN.equTime+tz*60)/1440
SUN.sunriseTime = SUN.solarNoon-SUN.HAsunrise*4/1440
SUN.sunsetTime = SUN.solarNoon+SUN.HAsunrise*4/1440
SUN.sunDuration = 8*SUN.HAsunrise
SUN.trueSolarTime = (SUN.time*1440+SUN.equTime+4*long-60*tz MOD 1440)
IF SUN.trueSolarTime/4<0 THEN
SUN.hourAngle = SUN.trueSolarTime/4+180
ELSE
SUN.hourAngle = SUN.trueSolarTime/4-180
ENDIF
SUN.solarZenithAng = DEG(ACOS(SIN(RAD(lat))*SIN(RAD(SUN.sunDec))+ _
COS(RAD(lat))*COS(RAD(SUN.sunDec))*COS(RAD(SUN.hourAngle))))
SUN.solarElevAng = 90-SUN.solarZenithAng

IF SUN.solarElevAng>85 THEN
SUN.atmosREfr = 0
ELSE
IF SUN.solarElevAng>5 THEN
SUN.atmosREfr = (58.1/TAN(RAD(SUN.solarElevAng))-0.07/ _
(TAN(RAD(SUN.solarElevAng))^3)+0.000086/(TAN(RAD(SUN.solarElevAng))^5))/3600
ELSE
IF SUN.solarElevAng>-0.575 THEN
SUN.atmosREfr = (1735+SUN.solarElevAng*(-518.2+SUN.solarElevAng* _
(103.4+SUN.solarElevAng*(-12.79+SUN.solarElevAng*0.711))))/3600
ELSE
SUN.atmosREfr = (-20.772/TAN(RAD(SUN.solarElevAng)))/3600
ENDIF
ENDIF
ENDIF
SUN.solarElevCorr = SUN.solarElevAng + SUN.atmosREfr
SUN.SolarAzAngle = (DEG(ACOS(((SIN(RAD(lat))* _
COS(RAD(SUN.solarZenithAng)))-SIN(RAD(SUN.sunDec)))/(COS(RAD(lat))* _
SIN(RAD(SUN.solarZenithAng)))))+180 MOD 360)
IF SUN.hourAngle < 0 THEN
SUN.SolarAzAngle = 360 - SUN.SolarAzAngle
ENDIF

END SUB

FUNCTION arrayOutput(arrayAz, arrayDec, arrayMax)
LOCAL az1, dec1, ac1, po

az1 = (SUN.SolarAzAngle-arrayAz+360 MOD 360)
dec1 = 90 - SUN.solarElevCorr
ac1 = DEG(ACOS(COS(RAD(dec1))*COS(RAD(arrayDec)) + SIN(RAD(dec1)) _
*SIN(RAD(arrayDec))*COS(RAD(az1))))
po = COS(RAD(ac1))*arrayMax

IF dec1 > 90 OR po < 0 THEN
po = 0
ENDIF
arrayOutput = po
END FUNCTION

sun_angle.zip

Jim

VK7JH
MMedit

lizby
Guru

Joined: 17/05/2016
Location: United States
Posts: 3399

Posted: 01:15pm 03 Aug 2025

Copy link to clipboard

Print this post

Thanks for that, Jim. Just an hour or so ago I asked perplexity.ai to compare fixed, daily, annual, and dual tracking at the equator and at 35 degrees latitude. Answer, at the equator, dual and annual tracking nearly worthless, daily tracking: 20-25% improvement over fixed. At 35°, a daily tracker adjusted manually adjusted 4 times a year improves over fixed by typically 25-30%. Only in places where you need to try to extract every watt is a dual-axis tracker likely to be called for.

PicoMite, Armmite F4, SensorKits, MMBasic Hardware, Games, etc. on fruitoftheshed

DaveJacko
Regular Member

Joined: 25/07/2019
Location: United Kingdom
Posts: 84

Posted: 09:01pm 03 Aug 2025

Copy link to clipboard

Print this post

Great work, interesting and informative..
Fortunately, in the North West of England, we don't need solar tracking.
Simply make sure the panel is the right way up

Try swapping 2 and 3 over

PilotPirx

Senior Member

Joined: 03/11/2020
Location: Germany
Posts: 101

Posted: 02:01pm 04 Aug 2025

Copy link to clipboard

Print this post

Hi Jim,

whats the meaning of these variables:
DIM array1Az = 90, array1Dec = 15 , array1Max = 2.5,
DIM array2Az = 270, array2Dec = 15, array2Max = 2.5,
Which variable contains the installation angle of the panels and their maximum output?

And how can I interpret the results of the day loop (the three values after the time, the last one is probably the sun angle)?

Peter

TassyJim

Guru

Joined: 07/08/2011
Location: Australia
Posts: 6297

Posted: 09:39pm 04 Aug 2025

Copy link to clipboard

Print this post

array1Az is the azimuth or direction the first array is pointing.
array1Dec is the declination or tilt of the first array in degrees
array1Max is the maximum power output of the first array.
The same for array2
I have two arrays both 2.5kW, one pointing east and one pointing west. Both on a 15 degree slope.

Quote PRINT thistime$,STR$(p1,2,3),STR$(p2,2,3),STR$(p1+p2,2,3), _
STR$(SUN.solarElevAng,3,1)

Time, power for first array, power for second array, total power, sun elevation.

I use a version of the code to produce this:

Jim
Edited 2025-08-05 07:51 by TassyJim

VK7JH
MMedit

Print this page

To reply to this topic, you need to log in.