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))