Solar Mirror Array Irradiance Calculator

--This assumes 100% reflective mirrors with M²=1 at all wavelengths. You'll likely need more to compensate for mirror inefficiencies.

local Temperature = 5777 --K
local MirrorRadius = 10 --m
local TargetDistance = 4.2e7 --m
local DesiredIrradiance = 1e8 --W/m²

local StartingWavelength = 1e-9 --m
local FinalWavelength = 20000e-9 --m
local Step = 1e-11 --m

local SolarIrradiance = 1367 --W/m²
--local MirrorEfficiency = 0.9 --unused
-- Yes, i know mirrors have different reflectances at different wavelengths, but i'm too lazy to actually look on material data for it.
--I'm stuck between using a silver(terrible at reflecting UV) or aluminum mirror(worse than silver on visible/IR light, but unlike silver it doesn't lose virtually all reflectance in UV).

local StefanBoltzmannConstant = 5.670367e-8 --W*m^-2*K^4
local BoltzmannConstant = 1.38064852e-23 --J/K
local PlanckConstant = 6.626070040e-34 --J/s
local SpeedOfLight =  299792458 --m/s

function PlanckIntegral(Wavelength, Temperature)
  --ported from http://www.spectralcalc.com/blackbody/appendixA.html
  local c1 =  (PlanckConstant*SpeedOfLight/BoltzmannConstant)
  local x =  c1  * (1/Wavelength) / Temperature
  local iterations = math.floor(math.min((20 + 100/x),1024))
  local sum=0
  for n=1,iterations do
    local dn=1/n
    sum= sum + math.exp(-n*x)*((x^3) + (3 * (x^2) + 6*(x+dn)*dn)*dn)*dn
  end
  local c2 = (2*PlanckConstant*(SpeedOfLight^2))
  return (c2*((Temperature/c1)^4)*sum*math.pi)/(StefanBoltzmannConstant*(Temperature^4))
end

local SpectrumCurve={}

print("Integrating blackbody spectrum...")
for n=StartingWavelength,FinalWavelength,Step do
  table.insert(SpectrumCurve,{n,PlanckIntegral(n,Temperature)})
end

print("Done. Calculating diffraction per wavelength...")
local TotalIrradiance = 0
for n=1,#SpectrumCurve-1 do
  local p = SpectrumCurve[n]
  local p2 = SpectrumCurve[n+1]
  local DiffractionAngle = 1.22*p[1]/MirrorRadius
  --print(DiffractionAngle)
  local BeamPower = (p2[2]-p[2])*SolarIrradiance*(math.pi*(MirrorRadius^2))
  --print(p2[2],p[2])
  local BeamIrradiance = BeamPower/(math.pi * (TargetDistance * math.tan(DiffractionAngle/2))^2)
  TotalIrradiance = TotalIrradiance+BeamIrradiance
  --print(string.format("Beam intensity between %3d and %3d nm is %e W", p[1]*1e9,p2[1]*1e9,BeamIrradiance))
end

print("Irradiance per mirror = "..TotalIrradiance.. "W/m²")
print(DesiredIrradiance/TotalIrradiance.."mirrors needed to achieve desired irradiance")