Here as a starting point is a quick cut-paste of some FORTRAN code for the zodiacal light that also includes a smooth version of the solar corona. An angular perturbation could be included for streamers. The 'vzod' result is the magnitudes per square arcsecond of the zodiacal light / corona, as seen from the vicinity of Earth. Other units can be derived such as spectral radiance and such.
-------------------------------------------------------------------------------------------------------------------------------
! Zodiacal Light S10 units 140 - 90 sin(beta) ecliptic latitude
elat = helioeclipticlat_a(ialt,jazi)
elon = helioeclipticlon_a(ialt,jazi)
if(elon .gt. 180.)elon = 360. - elon
s10zod_pole = 60. ! high ecliptic latitude
s10geg_ecliptic = 40. * cosd((elon-180.)/2.)**150.0
s10zod_ecliptic = 900. * cosd((elon)/2.0001)**7.7 + 140. + s10geg_ecliptic
! s10zod_ecliptic = 4000. * cosd((elon)/2.)**16.0 + 150. +
! s10geg_ecliptic
! Note the gegenschein of s10 = 40 can be added near
! helioecliptic 180.
eclp = 1.9 - 0.9*abs(elon/180.)
eclp = 1.0 ! control ecliptic sharpness
ecllatterm = abs(sind(elat))**eclp
s10zod = 10.**(log10(s10zod_pole) * ecllatterm + &
log10(s10zod_ecliptic) * (1.-ecllatterm))
! s10zod = (s10zod_pole) * ecllatterm + &
! (s10zod_ecliptic) * (1.-ecllatterm)
! s10zod = s10zod*10. ! debug for visibility
! F/K corona section (inner)
!
http://arxiv.org/pdf/0909.1722.pdf (far corona in Fig 1.)
!
https://ase.tufts.edu/cosmos/print_images.asp?id=28!
https://www.terrapub.co.jp/journals/EPS/pdf/5006_07/50060493.pdf! "Interplanetary Dust" edited by B. Grun
sqarcsec_per_sqdeg = 3600.**2
size_glow_sqdg = 0.2 ! sun/moon area
delta_mag = log10(size_glow_sqdg*sqarcsec_per_sqdeg)*2.5
distr = sqrt(elat**2 + elon**2)
distr2 = max(distr,0.40) ! account for pixel size
srad = distr2/0.25
spowerk = 5.7 + srad * 2.0
spowerf = 7.5 + log10(srad) * 2.0
! spowerf = 8.0 + log10(srad) * 2.0
! spowerf2 = 8.3 + srad * 0.16
smag_eff = -26.7 + 2.5 * min(spowerk,spowerf)
rmag_per_sqarcsec = smag_eff + delta_mag
elong = min(sqrt(elon**2 + elat**2),180.)
coeff = 3.5 - 2.5*cosd(elong/2)**180.
if(elong .gt. 0.)then
! pcoeff = 1.0 - elon**2/(elon**2 + elat**2)
! pcoeff = abs(elat) / sqrt(elon**2 + elat**2)
! The last exponent controls ecliptic sharpness
pcoeff = (abs(elat) / sqrt(elon**2 + elat**2))**1.3
else
pcoeff = 0.
endif
fexp = -2.5*(1.-pcoeff) + (-2.8*pcoeff)
f_wm2srum = coeff * srad**fexp
f_s10 = wm2srum_to_s10(f_wm2srum)
! f_nl = wm2srum_to_nl(f_wm2srum)
! Merge original zodiacal light with corona
s10zod2 = max(s10zod,f_s10)
sbu = s10zod2 * 1.25
vzod = s10_to_magsecsq(s10zod2)
glow_zod = log10(v_to_b(vzod))